Structure Formation in Modified Gravity Scenarios 
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We study the growth of structures in modified gravity models where the Poisson equation and the 
relationship between the two Newtonian potentials are modified by explicit functions of space and 
time. This parameterisation applies to the f(R) models and more generally to screened modified 
gravity models. We investigate the linear and weakly nonlinear regimes using the "standard" per- 
turbative approach and a resummation technique, while we use the spherical dynamics to go beyond 
low-order results. This allows us to estimate the matter density power spectrum and bispectrum 
from linear to highly nonlinear scales, the full probability distribution of the density contrast on 
weakly nonlinear scales, and the halo mass function. We analyse the impact of modifications of 
gravity on these quantities for a few realistic models. In particular, we find that the standard one- 
loop perturbative approach is not sufficiently accurate to probe these effects on the power spectrum 
and it is necessary to use resummation methods even on weakly nonlinear scales which provide the 
best observational window for modified gravity as relative deviations from General Relativity do 
not grow significantly on smaller scales where theoretical predictions become increasingly difficult. 

PACS numbers: 98.80.-k 



I. INTRODUCTION 



The discovery of the acceleration of the expansion of 
the Universe cannot be explained using General Relativ- 
ity and a matter content comprising only fluids with a 
positive equation of state. Seemingly, a new fluid with 
a negative equation of state, either a cosmological con- 
stant or dynamical dark energy, is required to generate 
the late time acceleration^]. Another plausible explana- 
tion could be that gravity itself is poorly understood on 
large scales and needs to be modified [2]. As General Rel- 
ativity (GR) is the unique Lorentz invariant low energy 
theory of spin two gravitons, any modification of grav- 
ity must include new degrees of freedom Q. Hence, in 
both the dark energy and the modified gravity contexts, 
new fields need to be included, the simplest ones being of 
course scalar fields. However, the presence of scalar fields 
is tightly constrained by fifth force and equivalence prin- 
ciple tests [J, This implies that the scalars leading to 
either dark energy or modified gravity must be screened 
in local and dense environments such as on earth or in 
the solar system B. Such models abound: chameleons @- 
III, dilatons [ilil - fl2l ]. GalileonsflUj. svmmetrons [l^ - [l6| and 
their generalisations [I?). In all these cases, the back- 
ground cosmology coincides with a ACold Dark Matter 
(ACDM) Universe. The only hope of observing non- 
trivial effects relies on the fact that perturbations in these 
models grow anomalously inside the Compton radius of 
the scalar field as first noticed in 0, [l8| . This anomalous 
growth can only be effective on intermediate scales. In- 
deed, on very large scales outside the Compton radius, 
normal gravity is retrieved while screening effects imply 
that GR is also recovered on small scales in very dense 
regions of the Universe Q- This opens up the possibility 
that relevant effects may be present at the mega parsec 



scale and that deviations from GR may be detectable by 
future galaxy surveys. 

In the following, we will concentrate on a formulation 
of the perturbation equations involving two Newtonian 
potentials and a time and scale dependent relationship 
between them. In terms of scalar field models, this cor- 
responds to the Jordan frame picture; the difference be- 
tween the two Newtonian potentials being due to the 
scalar field perturbation. In this picture, we choose to 
capture the modified gravity effects using a single func- 
tion e(k,a) whose interpretation in the Einstein frame 
is obvious: it measures the deviations of the geodesies 
under the influence of the scalar field. This function is 
universally characterised in terms of the mass and the 
coupling function of the scalar field. Here, we will con- 
sider it as defining the modified gravity models which we 
will study. 

Doing so, we neglect the nonlinear effects due to the 
presence of non-linear terms originating from the scalar 
field modifying gravity. As such we only modify the Eulcr 
equation by including the effects of a new scalar force. 
Hence, at this level of approximation, the models only 
differ from the GR treatment of ACDM perturbations by 
the inclusion of a time and scale dependent contribution 
to Newton's constant in the Euler equation. This simple 
modification of gravity is amenable to a quasi-linear and 
a fully non-linear treatment. 

The precision that future galaxy surveys will reach im- 
plies that simple linear perturbation theory is not accu- 
rate enough. One must include higher order effects and 
at one-loop order (i.e., next-to-leading order) we will find 
that the "standard" perturbative expansion is not suffi- 
ciently accurate to probe the modified gravity effects we 
investigate here. Therefore, we generalise a method de- 
rived using the saddle point of the generating functional 
of matter and velocity fluctuations. This resummation 
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scheme was already tested in the GR case and shown to 
be more accurate than the standard approach. 

To go beyond these low-order results we also study the 
dynamics of spherical perturbations floM^Tj ] . This can be 
exactly solved until shell crossing and it provides the full 
probability distribution of the matter density contrast on 
weakly nonlinear scales as well as the large-mass tail of 
the halo mass function. The latter can then be used to 
build a phenomenological halo model that also converges 
to the perturbative results on quasi-linear scales. This 
provides a simple estimate of the matter density power 
spectrum and bispectrum from linear to highly nonlinear 
scales, and a global picture of structure formation in such 
modificd-gravity scenarios. We discuss the relative devi- 
ations from GR of these various quantities as a function 
of scale. 

However, let us note that the analytical treatment of 
modified gravity developed here should only be taken as 
a first step, to indicate the type of effects one may ex- 
pect, because of our simplified parameterization of mod- 
ified gravity. First, more accurate modelizations would 
include some of the non-linearities due to the scalar po- 
tential at the one loop level[22j], which modify the Euler 
equation in an effective way. Second, the screening effects 
of the scalar field force in dense regions would mo dify 
the spherical collapse of an initial over density [13, l23j |. 
Here and as a first step, we will not consider these is- 
sues and treat the modification of gravity at the linear 
level in the scalar sector of the models. In (23 - [26| . this 
corresponds to the "no-chameleon" regime which should 
be seen as a non-screening case here in as much as we 
are neglecting the screening effects of modified gravity in 
dense regions. In the appendix, we compare our analytic 
treatment of the " no- chameleon" case with the simula- 
tions of [25} which shows a very convincing agreement. 
Of course, in future work, we intend to include one-loop 
corrections in the scalar sector as well as screening effects 
in the spherical collapse. Yet, it is useful to first develop 
the analytic formalism for the simpler parameterization 
studied in this paper. This will serve as a basis for more 
complex models that involve further ingredients (which 
are also more model-dependent, while the formalism de- 
veloped here can be applied to any function e(k, a) in the 
Euler equation). 

A similar approach was followed in [27j where f(R) and 
DGP models where considered. These cases were treated 
in the Jordan frame where the effect of modified grav- 
ity appears, for instance, in the difference between the 
two Newtonian potentials due to the anisotropic stress 
resulting from the presence of an extra scalar degree of 
freedom. In this work, the non-linear terms up to third 
order in the scalar dynamics were included, allowing one 
to study the onset of the screening mechanism at the 
perturbative level. Moreover, only the standard one loop 
contribution was taken into account in the quasi-linear 
regime and a fitting PPF formula was used to analyse 
fully non-linear scales. In the present work, the non- 
linearities in the scalar sector are not taken into account. 



On the other hand, we go beyond the standard one loop 
perturbative expansion and include a partial resumma- 
tion of perturbation theory. Moreover, the highly non- 
linear regime is studied using the spherical collapse and a 
halo model taking into account shell coupling due to the 
scale dependence of modified gravity. One of the advan- 
tages of our approach resides also in its versatility. Indeed 
we work in the Einstein frame where numerous models 
of modified gravity are defined [ljj]- Our treatment can 
be applied to chameleon and f(R) models and easily ex- 
tended to other models like dilatons and symmetrons. 
These extensions are being currently investigated. 

The paper is arranged as follows. In section UH we de- 
scribe the modified gravity models we will consider. We 
present the dynamical equations in the hydrodynamical 
approximation in section IIII[ and we study the pertur- 
bative regime in section IIV1 for the density power spec- 
trum and bispectrum. Next, we analyse the spherical col- 
lapse in the no-screening case in section [Vj This allows 
us to obtain the probability distribution of the density 
contrast on weakly nonlinear scales in section I VII and 
the halo mass function in section IVII1 Finally, we use 
these ingredients to build a phenomenological halo model 
in section IVIII1 which provides estimates of the power 
spectrum and bispectrum from linear to highly nonlinear 
scales. We conclude in section HXl 



II. MODIFIED GRAVITY 

A. The perturbed equations 

We consider models of modified gravity which can be 
defined by a change of the perturbation equations for 
Cold Dark Matter (CDM). The modifications are usually 
parameterised by two time and scale dependent functions 
7(fe, a) and (i(k, 0)12811. Other approaches have also been 
emphasized like in [29j . The 7 — fj, parameterisation does 
not follow directly from a Lagrangian formulation where 
causality is automatically taken into account. In the fol- 
lowing, we will use a restricted class of modified gravity 
models where the perturbed dynamics can be entirely 
specified by two time dependent functions only, m(a) 
and /3(a). These two functions enter as building blocks 
of a time and space dependent function e(fc,a). Finally, 
the knowledge of e(fc, a) defines 7(A), a) and /J,(k,a) com- 
pletely. The origin of this parameterisation springs from 
modified gravity models where a scalar field alters grav- 
ity on large scales and is screened in dense environments, 
leading to no modification of gravity in the solar system 
and in laboratory experiments. In turn, the dynamics of 
these models can be entirely reconstructed from the time 
evolution of the mass function to (a) of the scalar field, 
and its coupling to matter particles (3(a). This way of 
describing modifying gravity applies to chameleons and 
f(R) models, symmetrons and dilatons. Here, we will 
simply use the {m(a), /3(a)} parameterisation as a way of 
unambiguously defining modified gravity models at the 
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level of the perturbations. 

At the linear level, the perturbation equations of the 
CDM fluid follow from the conservation of matter 



(D 



where the density contrast is S = (p m — P m )/p m and 9 = 
d l Vi is the divergence of the velocity field. We denote by a 
prime the time derivative in conformal time t, with dr = 
At /a and a{t) is the scale factor. The Euler equation 
involves the Newtonian potential "J and reads in Fourier 
space as 



9' +H8 = fc 2 *, 



(2) 



where we denote Fourier-space quantities with a tilde. 
Here % = a' /a is the conformal expansion rate and we are 
using the Newtonian gauge with two distinct potentials 
<J and $, 



ds 2 



-a 2 (l + 2*)dr 2 + a 2 (l - 2$)dx 2 



(3) 



where x are comoving coordinates. The gravitational 
dynamics determine the evolution of 3> as 



- k 2 <b = Anu(k, a)Gp m 6/a, 



(4) 



which is a modification of the Poisson equation (p m is 
the mean comoving matter density and Q is Newton's 
constant). We also assume that there is a constitutional 
relation between the two potentials, 

* = 7(fc, a)l>, (5) 

implying that 

- fc 2 § = 47r ) u(fc, a)gp m S/a, (6) 

where 

jti(fc, a) = j(k, a)i/(k, a). (7) 
As a result, this implies that the density contrast obeys 



US' 



-H 2 p(k,a)5 = 0, 



(8) 



where Q m (a) is the matter density cosmological param- 
eter. The growth of structures depends on the choice of 
the function p,(k, a). We will define a large class of such 
models in the following section. 



B. Parameterised modified gravity 

The choice of function p,(k,a) seems to be unlimited. 
Here we focus on the simple choice 



and 



fi(k, a) = 1 + e(k, a) 
l + e(k,a) 



j{k,a) 



1 — e(k, a) ' 



(9) 



(10) 



where e measures the deviation from General Relativ- 
ity and is defined by two time dependent functions only, 
m(a) and /3(a) [13]. In modified gravity models with a 
screened scalar field in dense environments, m(a) is the 
mass of the scalar field at the cosmological background 
level. Similarly /3(a) is the coupling function between 
the scalar field and CDM particles. The space and time 
dependent function e(k, a) is expressed as 



e(fe, a) = 



2/3 2 (a) 



1 



m 2 (a)a 2 
P 



(ii) 



This parameterisation is valid for chameleons and f(R) 
models, symmetrons and dilatons[l7|. This implies in 
particular that 



p(k,a) = 



(1 



2/3 2 )fc 2 +m 2 a 2 



k 2 + m 2 a 2 



and 



7(fc,a) 



(1 + 2/3 2 )fc 2 + m 2 a 2 
(I - 2/3 2 )fc 2 +m 2 a 2 



(12) 



(13) 



This is an explicit parametrisation which shows that 
modified gravity effects only appear on scales such that 
k > amici), i.e. when scales are within the Compton 
wavelength of the scalar field. Outside the Compton 
wavelength, General Relativity is retrieved. These ex- 
pressions are valid in the Jordan frame where Newton's 
constant become time dependent too [I?} • For the models 
we consider here with m 3> H, such a time variation can 
be safely neglected in the Jordan frame. In the Einstein 
frame, the particle masses vary accordingly in a negligible 
manner. 

In the rest of this paper, we will only deal with one 
particular family of models defined by the coupling con- 
stant 
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V6 



(14) 



and the mass of the scalar field which is given by 

m(a) = m a- 3Cn+2)/2 , (15) 

where mo is a free scale which will be chosen to be close 
to 1 Mpc -1 and n > 0. In the matter dominated epoch, 
these models are equivalent to f(R) theories in the large 
curvature regime[17| where the f(R) correction to the 
Einstcin-Hilbcrt action reads 1 301 



f(R) w -167r£p A 



/go ^0 

n R" 



l+n 



(16) 



and p\ is the effective dark energy in the late time Uni- 
verse. In the recent past of the Universe, the mass of the 
large curvature models differs slightly from ([15]). see the 
appendix for more details. The mass mo is given by the 
useful relationship 



m 



Ho I ^mO + 4£!ao 
c 




i)I/hoI 



(17) 
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with c/Hq w 4 Gpc. Modifications of gravity must sat- 
isfy m c/ Hn> 10 3 to comply with a loosely screened 
Milky Way [31]. This also corresponds to |/,r | less than 
10~ 5 , the case \fno\ = 10~ 4 being marginal. When mo is 
too large, effects of modified gravity on large scale struc- 
ture occur on very non-linear scales. In the following, we 
will use values of mo ~ IMpc -1 which satisfy the loose 
screening bound for the Milky Way and imply interesting 
effects on large scale structure. 

We can also deduce now the two parametric functions 



and 



fj,(k, a) 



j(k, a) 



4 

3 in f. 



k 2 



4 k z 
3 rnr. 



2 fe 2 

3 rrif. 



where 



s = 3n + 4. 



(18) 



(19) 



(20) 



We will use the parameterisation of e[k, a) in the 
following when we give numerical examples. More 
precisely, we will consider the four cases (n, mo) = 
(0, 0.1), (0, 1), (1, 0.1), and (1,1), where m is given in 
units of Mpc -1 . This corresponds to the two scales 
mo =0.1 and 1 Mpc -1 and to the two exponents n = 
and 1. For these models we should have n > 0, see 
Eg. (1161) . and the choice n = for our numerical com- 
putations is only meant to exemplify the case of small 
n, that is s — > 4. The scales we consider are of the 
same order as the ones used so far in N-body simulations 
where \ f R0 \ = 10~ 4 , 10~ 5 , 10" 6 and n = 1. We will give a 
qualitative comparison with these numerical results, es- 
pecially we will briefly analyse the difference between the 
full numerical simulations, the no-chameleon case where 
the chameleon effects in dense region is neglected and our 
resummation method in the appendix. There we analyse 
the f(R) models where we take into account the late time 
effect of the cosmological constant on the mass function 
m(a). A more quantitative comparison is left for future 
work. 



III. PERTURB ATIVE DYNAMICS 

A. Hydrodynamical perturbations 

As explained in the previous section and in the intro- 
duction, we consider models where the continuity and 
the Euler equations are only modified by the non-trivial 
relationship between the two Newtonian potentials. For- 
mally, these equations have the same structure as in 
GR. When interpreted in terms of scalar field models, 
new non-linearities should appear in the Euler equation. 
However, the analysis of their role is left for future work. 



Then, the continuity and Euler equations read in Fourier 
space as 



g(k,r) + 0(k,r) 



dkidk 2 M k i +k 2 - k) 

xa(k 1 ,k 2 )0(k 1 ,r)J(k 2 ,r), (21) 



j£ (k, r) + H8(k, r) + ^H 2 [l + e(k, r)]S(k 7 r) = 

-J dkxdk 2 fe(ki +k 2 - k^k^kaJfltki.T^ka.-r), 

(22) 

which are the nonlinear generalizations of Eqs.([T]) and 
([2]), with the parameterization (|9]). The kernels a and (3 
are given by 



a(ki,k 2 ) = 



k\ 



(kx+kaj-ki |k 1+ k 2 | 2 (k r k 2 ) 
-,/3(ki,k 2 ) — 5p . 

In this paper we are mostly interested in the recent Uni- 
verse on large scales, hence we do not distinguish between 
the dark matter and the baryons that are treated as usual 
as a single collisionless fluid. These equations are only a 
first approximation of the dynamics of modified gravity 
on sub-horizon scales. Indeed, non-linearities in the po- 
tential and coupling function of the scalar field inducing 
the modification of gravity imply that the full dynamics 
should be described by the fluid equations for CDM par- 
ticles and the Klein-Gordon equation for the scalar field. 
Here we consider only the linear part of the scalar field 
dynamics which is tantamount to treating the scalar field 
as massive with a linear coupling to matter. When the 
mass of the scalar field is large enough m(a) 3> H, this 
allows one to integrate out the scalar dynamics and re- 
duce the equations of motion to the previous ones with 
a modified Newton constant. A priori, this procedure 
can be carried out to all orders taking into account the 
higher derivatives of the scalar field potential and cou- 
pling function at the minimum of the effective potential 
describing the background cosmology. Explicitly, this has 
been carried out to the one-loop level in the scalar field 
perturbation, resulting in an effective dynamics, once the 
scalar field effects have been integrated out, with a mod- 
ified /3(ki,k 2 ) [Hj]. The effect of this new contribution 
will be taken into account in a forthcoming publication. 

It is convenient to write the two fields S and 8 as a 
two-component vector tp [32[ ; which we define as 



tp2 



6 

-6/a 



(24) 



Because of the factor e(fc,r) in the Euler equation 
(|22l) the linear growing mode D + (k,t) depends on the 
wavenumber k. Therefore, instead of using D + as the 
time coordinate we use the logarithm of the scale factor, 



r]{t) = lna(i). 



(25) 
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This agrees with the standard choice used in most per- 
turbative studies for the simpler case of the Einstein- 
de-Sitter universe, where D + = a [32M35I ] . Then, the 
equations of motion (f2~T]) - (j2"2")) read as 



dip 

-Q^~i>i = I dkidk 2 5r>(ki+k 2 -k)a(ki,k 2 ) 



x ^ 2 (ki)-0i(k 2 ), 



(26) 



drj 



-Q m (l + e)ipx + 



1 



dc 



4'2 



2 v ,T \2 2 
dki dk 2 S D (ki + k 2 - k) (3 (ki , k 2 ) V 2 (ki ) V 2 (k 2 ) ,(27) 



where f2de(a) is the dark energy cosmological parameter 
and w the dark energy equation-of-state parameter. As 
in [33|, HH, H3| i this can be written in a more concise form 
as 

0(x, x') ■ i>{x') = K s (x; Xl ,x 2 ) ■ ^{ Xl )^{x 2 ), (28) 

where we have introduced the coordinate x = (k, r], i), 
i = 1, 2 is the discrete index of the two-component vector 
■ip, and repeated coordinates are integrated over. The 
matrix O reads as 



0(x,x') = S D (k-k r )6 D { v - v ') 



d_ 



-1 



(29) 



-ln m ( V )(l + e(k,r,)) 9_ + l-l w n d M 

and the symmetric vertex K s is 

K s (x;x 1 ,x 2 ) = Sr>(k. 1 +k 2 -k)S D (r)i-r))Sr)(r]2-'n) 
x 7 ^ ii2 (k 1 ,k 2 ), (30) 

with 

q(k 2 ,ki) a(ki,k 2 ) 
7i ; i,2(ki,k 2 ) = , 7 1;2il (ki,k 2 ) = , 



7 2 ; 2 , 2 (ki,k 2 ) = /3(ki,k 2 ), 



(31) 



and zero otherwise. 

The vertex K s does not depend on cosmology and it 
is not modified. Here modified gravity only affects the 
linear operator O through the term e(k, rf). In the case of 
a ACDM universe, that is, for e = 0, the matrix O and the 
linear growing mode D + (t) only depend on time. Then, 
it is possible to remove the explicit time-dependence of 
the equations of motion by using the time-coordinate 77 = 
ln£> + and making the approximation f2 m / f 2 — 1) where 
/ = dlnD + /dlna. This is a good approximation that 
is used in most perturbative works and it means that 
terms of order n in perturbation theory scale with time 
as D+{t) n [38j. Here we do not use this approximation 
because we consider the case where the linear growing 
mode and the matrix O also depend on wavenumber. 
This also means that in the ACDM limit, e — > 0, our 
approach is exact in the sense that it does not rely on 
the approximation fi m // 2 ~ 1. 



B. Linear regime 

1. Linear growing and decaying modes 

The linear regime corresponds to the linearization of 
the equations of motion (|2"51) or (|2"6")l - (|2T|) . We have al- 
ready discussed the linear equations in section Hi Al to in- 
troduce modificd-gravity effects. Here we present a more 
detailed analysis. The linear equations are O ■ i^l = or 



dip 



Ll 



drj 



i>L2 - 0, 



(32) 



^ - ^ m (l + e)i> L1 + Q - ^O de ) 4> L2 = 0, (33) 

where the subscript "L" denotes the linear solutions. 
Substituting Eg. (|3"2l into Eq. (|3"3")l yields a second-order 
equation for the linear modes D(rj), 

d 2 D (1 3 \ 3D 3 , 

v + (rHW^ M1+f)fl = a (34) 

As usual, we have a growing mode D + (rf) and a decaying 
mode D_(i]), and we define the initial conditions by the 
growing mode D + , so that in the linear regime we have: 



^i(k,?7) = <5xo(k) 



D+(k,ri) 



(35) 



In other words, we assume the decaying mode has had 
time to decrease to a negligible amplitude, which is the 
case in standard cosmologies. Then, the initial conditions 
are fully determined by the linear density field 5x,o(k). 

It is convenient to normalize the growing mode to the 
scale factor at early times. Indeed, we consider modificd- 
gravity models parameterized by a function e(k,a) such 
that e — > for a — > 0. Then, at early times we recover 
the Einstein-de Sitter universe (the dark energy compo- 
nent also becomes negligible) and we have the usual be- 
haviours: 

t^0: D+^a = e ri , D_ cx a~ 3/2 = e~ 3r ' /2 . (36) 

For numerical computations, it is convenient to introduce 
the reduced growing mode g+(k,7]) = D + (k,rf)/a. From 
Eq.® it obeys 



d 2 g+ /5 3 ^ \ dg+ 3 



d V 2 ' \2~2 WQde ) ~d^ + 2 K 1 - u; ) rj dc-^me].9+ =0 

(37) 



with the initial conditions 



7] -> -00 : g + -> 1, — > 0. 

or] 



(38) 



The linear growing mode can be easily computed from 
Eqs. (|3"Tl) - (f3"5|) . Although the linear decaying mode £)_ 




FIG. 1: Linear growing mode D+(k, t) normalized to the scale 
factor a(t) for four (n, mo) models. In each case we show 
the results for wavenumbers k — 1/iMpc -1 (lower curve) and 
k = 5/iMpc" 1 (upper curve), as a function of a(t). These two 
scales are in the non-linear regime and have only been chosen 
to exemplify the type of effects obtained in modified gravity. 



n = o, m„=o.i 



n — O, m n = 1 




n=l, m„ = 0.1 



n= 1 , m D = 1 

J I , I I L 



0.2 



0.4 0.6 
a 
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FIG. 2: Linear decaying mode D- (k, t) normalized to a(t)~ 3 ^ 2 
for four (n, mo) models. In each case we show the results 
for wavenumbers k = 1/iMpc -1 (upper curve) and 5/iMpc -1 
(lower curve), as a function of a(t). These two scales are in 
the non-linear regime and have only been chosen to exemplify 
the type of effects obtained in modified gravity. 



also obeys Eq. it is not convenient to use this for nu- 
merical computations (solving forward in time is unstable 
because of the contamination by the growing mode). It 
is better to use the Wronskian, 



W = D 4 



dD. 



dD, 



(39) 



drj drj 

which in our case is still independent of k and given by 



W{rj) 



-(1/2) /JW [l-3u>n dc (,,')] 



(40) 



FIG. 3: Linear growing mode D+(k, t) normalized to the scale 
factor a(t) for four (n, mo) models, at redshift z = up to 
non-linear scales. 
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FIG. 4: Linear decaying mode D- (k, t) normalized to a(t)~ 3 ^ 2 
for four (n, mo) models, at redshift z = up to non-linear 
scales. 



This normalization of W also defines the normalization 
of -D_, which reads 



J n D + {k,T]'Y 



(41) 



The integrals in Eqs. (|4"0")) and (|4T1) allow a fast computa- 
tion of D_(k, rj). 

We show in Figs.Q]and[5]the linear growing and decay- 
ing modes as a function of time (described by the scale 
factor a(t)). The deviation from the General Relativ- 
ity linear mode (which is almost identical to the lower 
curve in Fig. [1] and to the upper curve in Fig. [2]) in- 
creases for higher wavenumber. On these scales, the ef- 
fects of modified gravity grow as we span the parameters 
(n,m ) = (1,1), (0,1), (1,0.1), (0,0.1). Indeed, as seen 
from Eqs. (fTBl) - (f!?0"l) . deviations from GR appear at lower 
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k= 1 h Mpc - 




FIG. 5: Linear growth rate f(k,z) = dlnD+/d\na for 
wavenumber k = 1/iMpc -1 , for four (n, mo) models. 



k for small mass mo and at earlier time for smaller n. 
We can see that a positive e(k, a) in the Euler equation 
(f2"2"j) leads to a larger growing mode D + and a smaller 
decaying mode D_. This can be understood from the 
fact that a positive e can also be interpreted as a larger 
effective Newton constant in Eq.©. This implies a faster 
development of gravitational clustering and both linear 
modes evolve faster than in the ACDM cosmology. 

These behaviours can also be seen in Figs. [3] and 
[J] where we show the linear modes as a function of 
wavenumber at redshift z = 0. Although we plot our 
results up to k = lOO/iMpc -1 to allow a clear separa- 
tion between different curves, values beyond 1/iMpc^ 1 
do not describe the true quantitative difference between 
the models for observables such as the power spectrum 
because they are in the nonlinear regime, which is not 
described by these linear modes. In addition, on small 
scales new "screening" mechanisms, which are not de- 
scribed by the equations of motion (j2"Tj) -(f22" j) . take place 
and lead to a convergence to General Relativity and to 
the ACDM predictions. In agreement with the param- 
eterization (fTTj) . the linear modes deviate from the GR 
result at a wavenumber k ~ mo (in the plots the values 
of mo are given in units of 1 Mpc -1 ). At high k the de- 
viation is larger for smaller n (whence smaller s) because 
modifications of gravity have had more time to affect the 
dynamics, see Ea. ([l~8|) . 



2. Linear growth rate 

We plot in Fig. [5] the linear growth rate f(k, z) as a 
function of redshift, defined as usual by 



f(k,z) 



din D+{k,a) 
d In a 



(42) 



Both the linear growing mode D + and the linear growth 
rate / depend on wavenumber and to avoid overcrowd- 



ing the figure we only plot our results for k = 1/iMpc -1 
(which is in the mildly nonlinear regime at z — 0). The 
ACDM prediction could not be distinguished from the 
results obtained for (n, mo) — (1,1) and (0,1) (lower 
curves). In agreement with Fig. [I] the larger linear grow- 
ing modes D + obtained for (n, mo) = (1,0.1) and (0,0.1) 
lead to larger growth rates /. The deviation associated 
with the case (n, mo) = (1,0.1) would be difficult to 
detect with future surveys such as Euclid but the case 
(n,wp) = (0,0.1) should give a clear signal (see Fig. 2. 5 
in |1). 

3. Linear correlation and response functions 

From Eq. (j3"5)) the linear two-point correlation of the 
vector whence of the linear density and velocity 
fields, reads as 



C L {xi,x 2 ) 



D +1 D +2 D +1 D' +2 

5 D (k 1 +k 2 )P L0 (fci) 

D' +1 D +2 D' +1 D' +2 

3-D., 



(43) 
(44) 



where D + i = D + (ki,r]i) and D' + 



-(ki,Vi) 



In Sect. IIVB 21 we will consider a perturbative re- 
summation scheme that goes beyond standard one-loop 
perturbation theory. It involves the response function 
(or propagator) defined as the average of the functional 
derivative 



/ P#C!) \ 
R(Xl, X 2 ) = ( — — — - ) 



W> 2 )/<; =Q : 



(45) 



where £ is a "noise" added to the right hand side of 
Eq. (|28|) . Thus, R(xi,X2) measures the response of the 
system at time r\\ to an infinitesimal perturbation at an 
earlier time 772- It also describes the "propagation" of 
infinitesimal fluctuations. By causality, it satisfies 

m<V2- R(x 1 ,x 2 )=0, (46) 

and it obeys the initial condition 

rji^-vt'- R{xi,x%) -> <5o(ki -k 2 ) S iui2 . (47) 

In the linear regime, where the equation of motion 
reduces to O ■ ipL = 0, the response function obeys 

m>V2- O-R L = 0. 

Using the initial condition (|47|) . this gives 

©(?7i - ?/2)<Mki - k 2 ) 



(48) 



D' +2 D- 2 - D +2 D'_ 2 
D' +2 D-i-D'_ 2 D+i D- 2 D +X -D +2 D- X 

D' +2 D'_ x -D'_ 2 D' +l D^D'^-D^DU 



|(49) 



which involves both the linear growing and decaying 
modes D + and £>_. Here Q(r]i — r] 2 ) is the Heaviside 
function, which ensures causality. 
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IV. PERTURB ATIVE REGIME 

The equation of motion (|28[) is nonlinear and it has no 
explicit general solution. Therefore, it is usually solved 
by perturbative methods, which are sufficient on large 
scales and at early times where the density and veloc- 
ity fluctuations are small. Within our parameteriza- 
tion, modified gravity only changes the linear operator 
O of Eq. (l29l) . through the factor e(k,rj). Thus, we keep 
the same quadratic nonlinearity as in General Relativ- 
ity, with the same vertex K s of Eq. (f3T))) . Therefore, we 
can use the same perturbative schemes as in standard 
cosmologies. 

We first describe the standard perturbative approach 
in Sect. IIV Al and next a more accurate resummation 
scheme in Sect IIV B~2l Here we only go up to "one- 
loop order" : our standard perturbative prediction only 
includes the linear and one- loop (i.e., next-to- leading) 
terms, while our resummed prediction only adds a par- 
tial resummation of higher-order terms. 

We follow the approach described in detail in 40] (see 
also HI 111). 

A. Standard expansion 

Since the equation of motion (|28[) is quadratic in -0, it 
can be solved through a perturbative expansion in powers 
of the linear solution ipr,, as 

oo 

${x) = J2i> (n) ( x )> with i> (n) <* (^lT- (50) 

n=l 

Substituting this expansion into Eq. (I28[) gives the recur- 
sion 

n-l 

O-jW =K s (x;x 1 ,x 2 )-^2^(x 1 )^ n - e \x 2 ), (51) 

e=i 

which allows to compute terms of increasing order, start- 
ing with = tpL. One usually writes the expansion 
((50)1 in terms of the density and velocity fields, as [38ll4l1| 

OO p 

S(k,r;) = / dk i- k 

n<5o(ki + .. + k„ — k) 

n=l ^ 

x F^kx,..,^) Mki)..Mk„), (52) 

and 

OO p 

§(k, n) = Y< dk i-- k «Mki + .. + k„ - k) 

n=l ^ 

x E s n {\L X ,..,\L n -ri) 5 L0 (k 1 )J L0 (k„), (53) 

where Slo is the linear density field at some cho- 
sen time, as in Eq. (|35[) . The symmetrized kernels 
F% and are obtained from the recursion ([ST]) . In 



General Relativity the time-dependence of these ker- 
nels factorizes as F£ oc D"F*(ki, .., k„) and E^ cx 
— q(d In D + /dt)D+E%Cki , ..,k„) upon using the approx- 
imation f2 m // 2 ~ 1 [381 ] . In our case, where the linear 
growing mode D + (k, 77) depends on wavenumber, there 
is no such factorization and one must solve for the ker- 
nels F*(ki, .., k„; 77) and E„ (ki, .., k n ; 77) for each time 77 
of interest. 

Finally, from the expansion (|50l) one obtains the two- 
point correlation as 

C{ Xl ,x 2 ) = MuMxa)) (54) 

+^(2)^(2)) + ... (55) 

where we can use Wick's theorem to perform the average 
over the initial conditions ipLO- In particular, up to one- 
loop order the density power spectrum reads as 

P(k, 77) = P troc (fc, 77) + F llo °P(fc, 77), (56) 

where p tree ) associated with "tree diagrams", also corre- 
sponds to the linear power spectrum, 

P trec (M) = P L (k, V ) = D + (k, V ) 2 P L0 (k), (57) 

while p n °°P : associated with "one-loop" diagrams, is also 
given by 

P llo °P(k,r]) = P^(k,r])+P^(k,rj), (58) 
using the notations of with (see also [H, HI, E3 ) , 

P^(k,r,)=6P L0 (k) J dk' P La (k')F°(k',-k',k; V ), 

(59) 

P (C) (M)=2 J dk'P L0 (fc')^ (|k-k'|)P 2 s (k',k-k';7,) 2 . 

(60) 

B. Path-integral formulation 

1. General formulation 

The standard perturbative approach recalled in 
Sect. IIV Al computes the density power spectrum, and 
more generally many-body correlation functions, by first 
deriving an explicit expression for the nonlinear field ip 
in terms of the initial field ipL, as in Eqs. ([501) and ([52")) - 
(1531) . up to some order, and second taking the Gaussian 
average over the initial conditions, as in Eq. (|55[) . 

It is possible to work in the reverse order, by first tak- 
ing the average over the initial conditions and second 
writing an expansion in terms of the many-body correla- 
tions. A well-known procedure in the context of plasma 
physics and the study of the Vlasov equation is to use 
the BBGKY hierarchy, which gives a recursion between 
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successive correlation functions that may be truncated at 
some order [42| . A similar approach has also been used in 
35] to study the formation of large-scale structures in the 
single- flow perturbative regime, as in Eqs. ([2"Tj) - (l2"2"j) . As 
described in [33|, [111 Sjj, an alternative approach, also 
used in field theory and statistical physics [43|, EH, * s 
based on a path-integral formulation. There, it is shown 
that the statistical properties of the nonlinear field tp, 
which are fully defined by the equation of motion (|28|) 
and the Gaussian initial conditions (|35[) . can be obtained 
from the generating functional 

Z[j\ = (e^) = J T>4>T>\ (3-#- s &>%, (61) 

where A(x) is a Lagrange multiplier and the action S[i/j, A] 
reads as 

S[j>, A] = A • {O ■ $ - K s ■ fij>) ~ • A 7 • A (62) 

Here Aj is the two-point correlation of the initial con- 
ditions, taken at a time r\i. This matrix disappears in 
the final equations when we take the limit r/j — > — oo. 
Whereas moments of the field ip generate the many-body 
correlations of the density and velocity fields, such as the 
density power spectrum P(k), moments that involve the 
auxiliary field A generate the response functions [111 H|| . 
In particular, we have 

(A) = 0, (AA) = 0, $( Xl )\(x 2 )) = R{x x ,x 2 ). (63) 

As explained in [3^, S|| , the standard perturbative re- 
sults of Sect. IIV Al can be recovered from the generating 
functional (f6"Tj) . Indeed, one can see at once from Eq. (|5Tj) 
that the expansion (|50[) is also an expansion over powers 
of the vertex K s , with cx K^~ l and F£ cx if™ -1 . 
Therefore, the standard expansion in powers of 8lo for 
ip, which leads to the usual expansion in powers of Plo 
for averaged quantities, such as the density power spec- 
trum (|56p. is identical to an expansion in K s . Then, this 
expansion can be directly obtained from Eq. (|61l) by ex- 
panding in the cubic part A • K s ■ ipip of the action (|62p . 
This gives an alternative expression of the expansion (|55[) 
in terms of Feynman's diagrams [8 lj. 

2. Direct steepest-descent expansion 

One interest of the expression (151]) is that it can also 
serve as the basis of other approximation schemes. Here 
we focus on the "direct steepest-descent" method de- 
scribed in [33l |40j . which is compared with numerical 
simulations for the density power spectrum and bispec- 
trum in [§3, SH • In this approach, instead of expanding 
the cubic part of the action to write Eg. ([61]) as a series 
of Gaussian integrals, one expands around a saddle-point 
(which depends on 7) as in a semi-classical or "large- 
ly" expansion [46], [47| . This yields the Schwinger-Dyson 



equations 

O-C = S-C + n-i? T , (64) 
OR = 6 D + Z-R, (65) 

for the nonlinear two-point correlation C and response 
R, where £ and n are "self-energy" terms (there are two 
"correlations" , C and R, and two "self-energies" , £ and 
n, because there are two fields, the physical field tp and 
the auxiliary field A). 

These equations are exact and define £ and n. The 
"direct steepest-descent" or "large-N" expansion scheme 
corresponds to writing the self-energy terms £ and n 
as series in powers of the linear correlation Cl and re- 
sponse Rl- Then, the order of the approximation is set 
by the order of the truncation chosen for these expan- 
sions of £ and II. Because the truncation is made on 
£ and II, rather than on C and R, this automatically 
yields a partial resummation of higher-order terms (e.g., 
formally R would be given by the highly nonlinear expres- 
sion (0 — E)^ 1 whose expansion in Plo contains terms 
of all orders as soon as £ contains at least one power of 
Plo)- As described in [3^, [4(1 Sll , the result obtained for 
the correlation C at a given order (e.g., at one- loop order 
as in this paper) agrees with the result obtained by the 
standard perturbative expansion at the same order, and 
only differs by additional higher-order terms (which are 
only partially resummed). 

Then, this "direct steepest-descent" scheme gives at 
the one-loop order 

£ lloop (x,y) = ±K s {x;x u x 2 )K s (z;y,z 2 )R L {x u z) 

xC L (x 2 ,z 2 ) 1 (66) 

n llo °P(x,y) = 2K s (x;x 1 ,X2)K a (y;y u y 2 )C L (x 1 ,y 1 ) 
xC L (x 2 ,y 2 ). (67) 

This corresponds to a one-loop diagram [H, H3, S3] and 
at this order E cx Plo while II cx P^ . Substituting into 
Eqs. ([Ml) - (p5|) gives the nonlinear correlation complete up 
to order P£ , as in (|56")) , with the addition of a partial 
resummation of higher-order terms. Equation (|64l) can 
be solved as 

C(x u x 2 ) =Rx C L (vi) x R T + i? ■ II • R T , (68) 

where the first product does not contain any integration 
over time, and we take r/i — > — 00. Thus, to compute 
the density power spectrum up to one-loop order within 
the direct steepest-descent resummation, we first com- 
pute the linear correlation Cl and Rl, given by Eqs. (|44]) 
and (|49[) . This provides the self-energies E and II from 
Eqs. (p5|) and ([BT]) . Next, we compute R by solving the 
integro-differential equation (1651) and C from the explicit 
expression (f68|). 

The formalism used for the ACDM cosmology still ap- 
plies to our modelization of modified gravity. However, 
the numerical computation is somewhat heavier. Indeed, 
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as described in [H, in the ACDM case, the approx- 
imation f2 m // 2 ~ 1 allows us to explicitly factor the 
time-dependence of the linear correlation and response 
functions, and of the self-energies. Here this is no longer 
possible, because of the arbitrary function e(k,rj) in the 
linear operator ([29]) . This makes the numerical imple- 
mentation slightly more complex, as we can no longer 
use these factorizations to simplify the algorithms and 
we must keep track of the complex dependence on time 
and wavenumber of all linear modes and two-point func- 
tions. However, the method remains exactly the same, as 
described above, and it is still possible to devise efficient 
and reasonably fast numerical codes. 

3. Recovering the standard one-loop results 

Since we compute the self-energies £ and n for the 
one-loop steepest-descent scheme, we can also use them 
to recover the standard perturbative expansion instead 
of using the standard procedure recalled in Sect. IIV Al 
Indeed, the solution of Eq. (p5|) can be written as the 
expansion over powers of £, 

R = R L + R L -Y,-R (69) 
= R L + R L -Y<-R L + R L -Y 1 -R L -Y,-R L + ... (70) 

Therefore, up to order Plq we can write 

R = RW (71) 

with 

i? (0) = R L , R {1) = R L • £ lloop • R L . (72) 

Then, from (1681) the two-point correlation reads up to 
order P| as 

C = C (1) +C (2) , (73) 

with 

C« =R L xC L ( m )xRl = C L , (74) 

and 

= R W x C L ( m ) xRl + R L x C L { m ) x i?« T 
+R L • n lloop • R T L . (75) 

This expression is equivalent to Eqs. ([56| -([60 |) for the den- 
sity power spectrum [40j |. Therefore, since we have al- 
ready computed £ and II we can compute the standard 
one-loop power spectrum through Eqs. (j74|) - (j75)) . instead 
of using Eqs. (|59| -(|60 |) . This avoids explicitly computing 
the n— point kernels F£ of the standard expansion (|52[) . 

A similar procedure, based on the closure approxima- 
tion (34[, which is equivalent (at one-loop order) to the 
"2PI" effective action method of [33| , was used in [27| to 
obtain the standard perturbative predictions for several 
modified gravity models. However, while [27| included 



quadratic and cubic nonlinearities in the scalar field, as- 
sociated with the onset of the chameleon mechanism, in 
this paper we only consider modifications to the Poisson 
equation at the linear level. On the other hand, within 
our simpler formulation of modified gravity we go be- 
yond the standard perturbative approach by computing 
the "steepest-descent" resummation presented in the pre- 
vious section. 



4- Alternative resummations 

Finally, the path-integral (|6 X |) can also lead to alterna- 
tive resummation schemes, such as the "1PI" and "2PI" 
effective action methods described in [46] . The 2PI effec- 
tive action still leads to the Schwinger-Dyson equations 
(|64l) - (|65[) but the self-energy terms are given in terms of 
the nonlinear two-point functions C and R, instead of the 
expansion over Cj, and Rl used in the direct steepest- 
descent scheme. At one-loop order, this amounts to re- 
placing C L and R L by C and R in Eqs.{66l)-(!67]). How- 
ever, already for the ACDM case this makes the compu- 
tation more complex since Eg s. ([Ml) - (155)) become coupled 
nonlinear equations over C and R [33t HH. Then, one 
needs to solve for the four quantities C, R, S, and n by 
simultaneously moving forward with time. This numer- 
ical computation was performed in [33[ and it appeared 
that it did not provide a significant improvement over the 
simpler direct steepest-descent scheme (although a more 
precise comparison with numerical simulations may re- 
main of interest). Therefore, we do not investigate this 
scheme further. 

The direct steepest-descent method of Sect. IIV B 21 is 
not necessarily the most accurate resummation scheme. 
In particular, it yields a response function that does not 
decay at high k or late times, but shows increasingly fast 
oscillations with an amplitude that follows the linear re- 
sponse function. This is not realistic, since one expects 
a Gaussian-like decay for Eulerian response functions, as 
can be seen from theoretical arguments and numerical 
simulations [H, Hfl li^ - IHoj . However, the fast oscilla- 
tions still provide an effective damping in a weak sense 
(that is when the response function is integrated over). 
Alternative resummation schemes have also been studied 
in the literature, such as the "renormalized perturbation 
theory" [H, HH and several related approaches [5ll - [53l j , 
which rely on a response function that interpolates be- 
tween its low-fc standard perturbative expression and a 
resummed high-fc limit, methods based on path-integral 
formulations [54|, on closures of the hierarchies satisfied 
by the correlation functions [H], |35[ , or on Lagrangian- 
space formulations [55j . 

The reason why we consider the direct steepest-descent 
method here is that it provides a simple and efficient 
method, which has already been shown to be reasonably 
accurate for ACDM cosmology [3?], HE)]. An advantage 
with respect to some alternative approaches, which can 
show similar levels of accuracy, is that it is fully system- 
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atic and contains no free parameter or interpolation pro- 
cedure. Therefore, the generalization from the ACDM 
cosmology to modified-gravity scenarios is straightfor- 
ward, as described in Sect. IIVB 21 and we can expect 
a similar accuracy. 

C. Bispectrum 

Because the gravitational dynamics is nonlinear, the 
density field becomes increasingly non-Gaussian in the 
course of time. The most popular measure of these 
non-Gaussianities, which can be used to break degen- 
eracies between cosmological parameters or to constrain 
primordial non-Gaussianities, is the three-point correla- 
tion function (56[. In Fourier space this is the so-called 
bispectrum, 

(5(ki)5(k 2 )*(ka)) = fe(ki+k 2 +ka)B(*i,fc2,*3). (76) 

This can be computed by the standard perturbative ap- 
proach (38|. Substituting the expansion (I5"2"j) yields the 
standard tree-order result 

fl^fci.fc.fc) = 2F°(k 2 ,k 3 ; V )P L0 (k 2 )P LQ (k 3 ) +2 eye. 

where "2 eye." stands for two terms obtained by circular 
permutations over {k\, k 2 , k 3 }. 

Within the path-integral formalism of Sect. IIVB 11 ex- 
panding Eq. (j61l) in powers of K Sl that is, in the cubic 
part of the action S, yields for the three-point correla- 
tion at tree-order [40| 

C* rco = Rl ■ K s ■ ClCl + 5 perm. (78) 

This gives for the equal-time density bispectrum: 

i? trcc (fci,fc 2 ,fc 3 ;?7) = 2 f d v ' Y, Rmafawrf) 

J —00 -/ -/ -/ 

x T 4. i2)i3 (k 2 ,k 3 ) + 2cyc. (79) 

which is again equivalent to Eq. (177[) . In practice, instead 
of Eq. ([77j) we use Eq. (j?9l to compute the standard tree- 
order bispectrum. The effects of the modified-gravity 
function e(fc, a) are included through the linear correla- 
tion and response Cl and Rl, which depend on the mod- 
ified linear modes D + (k,a) and D_(k,a) as described in 
Sect. IIII1 As in Sect. IIV B 31 this allows us to obtain the 
"standard" perturbative predictions without computing 
the kernels F£ of Eq.(j52]). 

At one-loop order the expressions involve more terms. 
They can be found in [4(| (for the ACDM cosmology) for 
the standard approach as in ([77)1 . the equivalent path- 
integral formulation as in (|79[) . and the direct steepest- 
descent method used in Sect. IIVB 21 for the power spec- 
trum. Contrary to the power spectrum, a detailed com- 
parison with numerical simulations [45| shows that at 



one-loop order the steepest-descent resummation for the 
bispectrum is not more accurate than the standard re- 
sult. Therefore, we do not investigate this resummation 
for the bispectrum here. 

Because the linear modes depend on wavenumber, 
computing the one-loop order terms is significantly more 
difficult than in the ACDM case, even within standard 
perturbation theory. Using the scalings B treo cx Z>iPj; 
and B lloop cx D+Pl , we consider the following approxi- 
mation: 

(Dtrec \ 3/2 
ACDM / 

Thus, we simply rescale the one-loop correction ob- 
tained in the ACDM scenario by the prefactor 

(B tICC /BXcb M ) 3/2 - This would be exact if the ratio 01 
the linear modes were constant. We choose this prefac- 
tor, rather than (D + (k)/D +i ACDM(k)) 6 , because it in- 
cludes an integration over the past history and over the 
appropriate range of wavenumbers of the linear modes. 
This should be sufficient for our purpose, which is simply 
to estimate the magnitude of these one-loop corrections. 



D. Numerical results 

1. Set up 

For our numerical computations, we adopt in 
this paper a flat ACDM reference model with 
cosmological parameters (^ mj h, cr 8 , n s ) = 
(0.279,0.046035,0.701,0.817,0.96), which is consis- 
tent with WMAP 5-year observations [Hj]. We use a 
publicly available code, CAMB [59j, to compute the linear 
power spectrum including baryon acoustic oscillations. 
This is the same cosmology as used in [§3, 45] , which al- 
lows a clear comparison with their ACDM results. Then, 
the four models that we consider in this paper, defined 
by the parameters (n, mo) = (1, 0.1), (1, 1), (0, 0.1), and 
(0, 1), as described in Sect. Ill Bl arc defined by the same 
initial conditions as this reference ACDM model. This 
means that they all coincide at early times and on large 
scales, because e(k,a) — > for a — > or k — > 0, but 
their linear variance as today on scale 8/i _1 Mpc slightly 
differs. 

For later use, let us note <5l(a) (x, 77) the linear density 
field within the reference ACDM cosmology, 

S L (A)(Kv) = D +(A) (r 1 )~S L0 (k), (81) 

where D+(A) is the ACDM linear growing mode, which 
does not depend on wavenumber. Then, the actual linear 
density field can be written in terms of this reference 
ACDM linear field as 

~S L {k, V )= ° +{k 'f: U z(A)(M). (82) 
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FIG. 6: Ratio of the power spectrum P(k) to a smooth ACDM linear power spectrum PLs(k) without baryonic oscillations, 
from [53]. We show our results for three models with (n, mo) = (1,0.1) (middle red lines), (0,0.1) (upper black lines), and 
(0,1) (lower blue lines). In each case, we plot both the linear power (dashed line) and our nonlinear result (solid line) from 
Eg . (]124[) . which is based on Eg . (|68[) . For comparison, we also plot the standard 1-loop result from Eg. (|73p for the case (0, 1) 
(upper blue dotted line). 
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FIG. 7: Ratio of the equilateral bispectrum, B cq (k) = B(k, fc, fc) , to the product 3Pi s (fc) 2 , where Pha{k) is a smooth ACDM 
linear power spectrum without baryonic oscillations, from [57]]. As in Fig. [5] we show our results for three models with 
(n, mo) = (1,0.1) (middle red lines), (0,0.1) (upper black lines), and (0, 1) (lower blue lines). In each case, we plot the tree- 
level bispectrum (dashed line) from Eq. (|79[l . the 1-loop bispectrum (dash-dotted line) from Eq. (|80[l . and our nonlinear result 
(solid line) from Eq.(fl28}. 



This is merely a re- writing of the initial conditions, which 
we choose to express at any time r\ through the reference 
ACDM growing mode. 



2. Power spectrum 

We show our results for the matter density power spec- 
trum P(k) on BAO (baryon acoustic oscillations [601 ]) 
scales in Fig. [5] To clearly distinguish the different curves 
and the baryon acoustic oscillations we normalize P{k) 
by a smooth ACDM linear power spectrum PL S (k) with- 
out baryon oscillations, from (5?]]. Our nonlinear pre- 
diction includes both the perturbative "two-halo" part 
P2H.(k), based on the steepest-descent resummation (1551) . 



and the nonperturbative "one- halo" part Pm(k), as de- 
scribed in Sect. I Villi and Eq. ([124[) below. However, on 
these scales the power spectrum is dominated by the per- 
turbative contributions and the full nonlinear result is 
very close to the resummed perturbative part (|68[) . 

As explained above, all our results converge at low k 
to the same reference ACDM power, P^(k), because 
of our common choice of initial conditions. Moreover, 
on the scales shown in Fig. [SI this ACDM power spec- 
trum cannot be distinguished from the (n = 0,mo = 1) 
result, where the effects of modified gravity are the weak- 
est amongst the models that we consider here. As in the 
ACDM cosmology, the nonlinear evolution amplifies the 
power spectrum but erases most of the oscillations. The 
difference between the various modified gravity models 
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and General Relativity is rather small and it is not am- 
plified by the nonlinear evolution. We clearly see that to 
probe these deviations it is necessary to go beyond lin- 
ear theory and to include at least one-loop corrections. 
Moreover, the comparison with the upper dotted curve, 
which shows the standard one-loop result for the case 
(n = 0,mo = 1) (which cannot be distinguished from 
GR), shows that these modified- gravity effects are at the 
order of or smaller than the accuracy of the standard 
one-loop prediction. This means that to probe modified 
gravity on these scales it is necessary to use more ac- 
curate analytical formalisms, such as the resummation 
scheme described in Sect. IIVB 21 and used in this paper, 
or to include higher-order corrections within the standard 
perturbative approach (but this latter option may not be 
very efficient because the standard perturbative expan- 
sion does not converge very well). This provides another 
motivation for the development of efficient perturbative 
schemes, which re-sum high-order contributions. 



3. Bispectrurn 

We show our results for the matter density bispectrurn 
on BAO scales in Fig. [7J Here we only consider equi- 
lateral configurations, B cq (k) — B{k,k,k), and we nor- 
malize the bispectrurn by 3Pj, s (fc) 2 . Because PLs(k) is 
not the actual power spectrum but a smooth ACDM lin- 
ear power spectrum without baryon acoustic oscillations, 
this ratio is not identical to the usual "reduced bispec- 
trurn" Q eq = B eq / (3P 2 ). However, this allows us to 
clearly distinguish the baryon acoustic oscillations of the 
tree-level bispectrurn (|77| - ([75|) . Again, on these scales 
the ACDM bispectrurn cannot be distinguished from the 
(ri = 0,rno = 1) result. 

As for the power spectrum shown in Fig. [51 the nonlin- 
ear evolution amplifies the bispectrurn but erases most of 
the oscillations. The difference between the various mod- 
els and GR is again rather small and it is necessary to 
go beyond the tree-level prediction. Unfortunately, the 
comparison between our approximate one-loop prediction 
and our full nonlinear model, which includes the non- 
perturbative "two-halo" and "one-halo" contributions as 
described in Sect. I Villi below, suggests that one-loop 
terms are not sufficient to obtain reliable measures of 
such modified-gravity effects and that nonperturbative 
contributions cannot be neglected. Since the theoreti- 
cal accuracy of such nonperturbative terms is lower than 
the one of perturbative terms (which can be computed in 
a systematic and rigorous fashion), this means that the 
bispectrurn is not a very efficient probe of these modified- 
gravity models (unless one can run dedicated N-body 
simulations for each modified-gravity scenario). Thus, 
the power spectrum studied in Sect. II VP 21 should pro- 
vide a better tool, as the accuracy of its theoretical pre- 
dictions is better controlled. 



V. SPHERICAL COLLAPSE 
A. General case 

To go beyond low-order perturbation theory, the main 
analytical tool that can provide exact nonlinear results 
is the study of the spherical collapse. This allows an ex- 
plicit computation of the nonlinear dynamics (restricted 
to spherical symmetry) that can also serve as a basis to 
evaluate several quantities of cosmological interest, such 
as the halo mass functions and the probability distribu- 
tions of the density contrast. We describe in this section 
the equations that govern the spherical dynamics and 
give a simple approximation for typical fluctuations. 

Following the usual approach for ACDM or 
quintessence cosmologies [6l|, [62j, the physical ra- 
dius r(t), which contains a constant mass M until 
shell-crossing, evolves as 



dr 



1 a* 

a dx ' 



with * = $ 



N 



(83) 



where ^ is the total potential seen by massive particles. 
Here we note with a dot derivatives with respect to time 
i, physical coordinates by r and comoving coordinates by 
x. Within our framework, defined by Eqs. (f!?l"j) - (|2"2"1) . the 
potential VE' contains two parts, the usual Newtonian po- 
tential <I>n = $n , associated with General Relativity, and 
the effective component \I> £ , associated with the modifi- 
cation of gravity. 

In physical coordinates, we have 



V 2 $ N 



4^(p^-) + (l + 3Hpl hys ' ) ), (84) 



where we note with a superscript "(phys.)" densities in 
physical coordinates and we again assumed an uniform 
dark energy component. Using Gauss' theorem, this 
yields the usual part (r)^ of the acceleration of the shell 
at radius r roll. 16211 . 



(On = - 



47r£ 



pL PhyS °«r-) + (l + 3^ p e hys ' ) l ) (85) 



where pm hys '''(< r) is the mean physical density within 
radius r, 



pL phys) (<0 



3M 
4 7rr 3 



(86) 



In comoving coordinates (with the background Hubble 
flow), the effective component VE^ only depends on the 
matter density fluctuations, Spm — p m — p m , through 

§ e = e(fc,*)5$N with V 2 ((5$ N ) = 4ng5p m /a, (87) 
whence 



* e (M) = -^^e(MWM)- 
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This is a linear approximation in the spherical collapse 
dynamics which is only valid as long as the screening 
effects of modified gravity are not taken into account. 
When the screening effects appear, the scalar force lead- 
ing to the extra contribution in Newton's equation is 
highly suppressed and the spherical over density collapses 
like in GR. These effects can be modeled out in the top- 
hat approximation like in [20| or using the exclusion set 
theory (23[. Taking into account these effects is left for 
future work. 

Going back to configuration space, this yields the ad- 
ditional part (f) e due to this "fifth force", 



the previous sections, Eq. (|93]) reads as 



-l!!ii r p(phy S .) f dk4irk 2 e(k)6(k)W(kx), 
3 Jo 

(89) 

where the integral is written in terms of comoving quan- 
tities and x ~ r I a. Here we introduced the Fourier trans- 
form of the 3D top-hat of radius x and volume V, 



W{kx) 



dx' iW sin(fcir) — kxcos(kx) 

- e = 3 V) 3 • (90) 



If e does not depend on wavenumber we can check that 
Eq.® gives 



4tt£ 

= — e r 

3 



Jphys. 



(<r) 



^(phys.) 



(91) 
(92) 



In agreement with Eq.(|85j). an uniform e(t) gives rise to 
a fifth force that is proportional to the Newtonian grav- 
itational force where we subtract the background part 
(associated with the mean density of the universe). 

Collecting Eqs.(|55|) and (|59")l we obtain the equation of 
motion 



4tt£ 



P^ s:) {<r) + {l + 'iw)pT yS - ] 



+p(P hys ) / dk4irk 2 e(k)5(k) W{kx) 



(93) 



As in |6lLl62j |. it is convenient to introduce the normalized 
radius y{t) defined as 



r(t) 

y(t) = -£f with q 
a(t)q 



3M 

^Pm 



1/3 



y(t = 0) = l. 



(94) 

Thus, q is the Lagrangian comoving coordinate of the 
shell r(t), that is, the comoving radius that would enclose 
the same mass M in a uniform universe with the same 
cosmology. This also implies 



pL phys) (< r) 

-^(Phys) 



y~ 3 , S r = 5(< r) = y- 3 - 1. 



(95) 



d^y 
drj 2 



1 3 \ dy 

2 " 2 Wndc ) ^ 



"-^-y / dk4Trk 2 e(k)5{k)W(kx). (96) 
2 Jo 



Choosing again rj = In a(t) as the time coordinate, as in 



The left hand side is the usual result in ACDM cosmology 
[6l],|62[ and the right hand side is the new term associated 
with the "fifth force" . If e does not depend on wavenum- 
ber, the integral reduces to e(a)5(< x) — e(y~ 3 — 1), as in 
the usual third term of the left hand side. Then, the mo- 
tion of each mass shell, described by y(M,r]) or r(M,r]), 
is independent of the other shells before shell crossing. If 
e(k, a) depends on wavenumber, the integral does not re- 
duce to a simple function of y at the same mass scale and 
it explicitly depends on the whole density profile, S(x) or 
S(k) in Fourier space, of the matter perturbation. Then, 
the dynamics of all mass shells are coupled at all times, 
even before shell crossing, and we must solve for the evo- 
lution of the full density profile with time, y(M, 77), as a 
function of M and rj. 

In previous works [U [24|], the spherical collapse dy- 
namics was often approximated through an effective 
rescaling of Newton's constant (this corresponds to a 
function e(a) that does not depend on k). This allows 
one to recover the usual form of the equations of motion 
where all shells are decoupled before shell crossing. By 
varying this effective Newton constant [HI, or making it 
a dynamical variable that depends on the environment 
ETil . one may capture screening effects. Here we do not 
include such screening effects but Eq. (|96]) takes into ac- 
count the dependence of the dynamics on the density pro- 
file. This allows us to include the effects associated with 
the dependence on wavenumber of e(fc,a). As we will 
check in Fig. |8] below, this already yields a dependence 
on mass of the linear density threshold 5 C (M) associated 
with halo formation. 

Thus, the modified-gravity term makes the equation 
of motion significantly more complex, because it is no 
longer local and it turns the usual ordinary differential 
equation into a partial integro-differential equation. 



B. Approximation for typical profiles 

Let us assume we are interested in the dynamics of 
a single mass shell M. Then, we wish to obtain from 
Eq. dM)) a closed approximate equation for yM(rj) = 
y(M, ij), which does not involve the other shells M'. The 
simplest method is to use an ansatz for the density pro- 
file S(x,r]), or S(k 7 r]), that is parameterized by yM(v)- 
This will simplify numerical computations because it will 
transform Eq. (I96p into a single ordinary differential equa- 
tion. Then, let us recall that the mean conditional profile 
of the linear density contrast 5l (x) , under the constraint 
that the mean density contrast within a comoving radius 
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R is equal to 6 lb., reads as 

fe(*) = ^r / ^<Wx,xO, (97) 

where Cs L s L is the matter density linear correlation, 
C , j i { i (xi,x 2 ) = (<5l(xi)5l(x 2 )) 

= fdHrfP4fe)?f^l, (98) 
7o fc|x 2 -xi| 

and a R is the variance of the linear density contrast at 
scale R, 

/>oo 

Vr = (S 2 lr)= AkATrk 2 P L {k)W{kRf. (99) 
Jo 

This only relies on the assumption that the linear density 
field is Gaussian. Then, we consider the approximation 
where the density profile used in Eq. (|96p is set to 

5(x) = %-/ p-C Sx . SL M, (100) 
a i M JVm vm 

which reads in Fourier space as 

8{k) = p dk) W{kx M ), (101) 

where we used S XM = yjf — 1. Here xm{v) — 
r M(v)/ a ( r l) = UMiVjQM is the comoving radius of the 
shell M and it follows its spherical dynamics. Substitut- 
ing the ansatz (|101[) into Eq. (|96| gives the equation of 
motion 

d 2 2/M (1 3 \ dy M O m , _ 3 , 

-dtfT + 1 2 " 2 ™" dc J + T {vm ~ 1} yM 

x (l + / dfc47rfc 2 e(fc)P L (fc)iy(fca; M ) 2 ) = 0. 

V a t M Jo J 

(102) 

The equation (|102j) is exact if e does not depend on 
wavenumber, in which case the parenthesis is equal to 
(1 + e(r))) and we recover the behaviour of Eq. ([92p . It 
is also valid at order one over Sl and e when the ini- 
tial perturbation has the linear profile (|97p at early time. 
Thus, it agrees with the typical profile (|97|) . under the 
constraint &lx m aTi mass shell M, in the linear regime, 
at zeroth order over e. It is no longer exact at higher 
orders over Sl because the nonlinear dynamics changes 
the shape of the density profile in a complex fashion. It 
is not valid at order e, even in the linear regime, because 
the mean profile (|97[) is not a solution of the linear dy- 
namics, as the linear growing mode D + (k, a) depends on 
wavenumber. In our case, where e <C 1, this is a negligi- 
ble effect and we would actually obtain similar results by 
using in Eqs. (|100[) and (|101j) the reference ACDM linear 
correlation Cs L s L (A) an d power Pl(A)- 
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FIG. 8: Reference linear density contrast 5 c (a) = J'q 1 (200) 
associated with a nonlinear density threshold of 200 at red- 
shift 2 = 0. We show our results as a function of the halo mass 
M for four (n, mo) models, for typical initial profiles of the 
form (1104[) . In each case, the upper curve is the approximate 
result from Eq.fTD5]) and the lower curve the exact result from 
Eq. (|96j) . 

C. Spherical-collapse mapping 

In the linear regime we can check that Eq. (l96[) agrees 
with Eq. (|34[) for the linear growing mode. Indeed, using 
Ul = 1 — SLq/3, SLq — J v dx(5z,(x)/V r , and x — q at lowest 
order, Eq.tjMf becomes at linear order: 

J V vJ dke {^ k ) + U"2 WQde JW (k) 

_^(l + e (A))^(k)j=0. (103) 

This agrees with Eq. ([M|) and we recover the linear solu- 
tion <5 L (k, i]) = D + (k, T])S L0 (k). 

At linear order, the ansatz (|10ip reads in Fourier-space 

as H ft) = ( S Lq M /crq M )PL(k)W(kq M )- Substituting into 
Eq. (|103p remains exact if the profile of the perturbation 
is given by Eq. (|100p (or for the shell M, whatever the 
initial profile, if e does not depend on wavenumber). 

We now consider the spherical dynamics of typical ini- 
tial perturbations, of the form (|97p at early times, which 
we write as 

2 

W(A)=<WA)T^ (104) 

^(A) 

for the mean initial density contrast within arbitrary ra- 
dius q'. Here, as explained in Sect. IIVD 11 we choose 
to write the initial conditions in terms of the reference 
ACDM linear field, which is simply an "update" at arbi- 
trary time r\ of the initial field Slo given at a fixed time. 
This is more convenient than using the actual linear field 
Sl , which depends on the modified-gravity growing mode 
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77) and mixes dependences on the initial conditions 
and on the modified gravity parameters. In this fash- 
ion, Eq. (|104j) describes the same initial condition for all 

our models. Here a 2 is the cross-correlation of the 

91,92 (A) 

smoothed reference linear density contrast at scales q\ 
and q 2 , 

a quq 2 W = (<W(A)<5z, fc (A)) 

= / dk4TTk 2 P L(A) (k)W(k qi )W(kq 2 ) 7 (lQ5) 
Jo 

and Vq/fA = c 2 q ^ ■ For each mass scale q, with M — 
(47r/3)p m g 3 , and initial amplitude 5 Lq ^, which define 
the initial condition (|104j) . we can solve the spherical dy- 
namics or the approximate dynamics (|102l) . For the 
"exact" dynamics (IM1) we consider for simplicity that in- 
ner shells that have already collapsed to the center of the 
halo remain at the center. (After shell crossing we should 
modify Eq. (I96p to take into account the change with time 
of the mass enclosed within a given shell. However, we do 
not consider this effect because radial orbits suffer from a 
strong instability, which diverges at the time of collapse 
to the center [6J] , and after that time one should include 
transverse motions that lead to virialization.) As long as 
shell crossing is restricted to inner shells, within the mass 
scale M of interest, this is not a very serious problem be- 
cause the dynamics is mostly sensitive to the total mass 
enclosed within a given radius (as in the usual Newto- 
nian case or for e that does not depend on wavenumber) 
or to the local slope of the density profile (for the \ow-k 
behaviour e(k) oc k 2 ). 

At a given mass scale q and time rj, this defines a map- 
ping, <5l 9 (a) h- »• 6 X = -^(^^(A)); from the reference linear 
density contrast SLq(A) to the nonlinear density contrast 
8 X . Here x is again the Eulerian comoving radius of the 
shell M, with x = r/a = yq as in (|94|) . 

If e does not depend on wavenumber, this mapping 
does not depend on the scale q nor on the shape of the 
initial profile. If e depends on wavenumber, this mapping 
depends both on the mass scale q (whence the subscript q 
in !Fq) and on the initial shape of the profile (which is why 
we had to choose a specific case, such as the typical shape 
fll04[) ) . This implies that if we choose for instance a given 
nonlinear density threshold, such as 200, to define halos, 
the associated linear density contrast 8 C ^ = J r 9 " 1 (200) 
depends on the mass of the halo (through the scale q) . 

We show our results for this linear density threshold 
J'~ 1 (200) at redshift z — in Fig. [SI For each model we 
plot both the exact result from Eq. flMl) and the approx- 
imate result from Eq. (|102[) . We clearly see the mass de- 
pendence associated with the modification of gravity. For 
positive e gravitational clustering is more efficient and a 
lower value of <5l(a) is required to reach the nonlinear 
density contrast 6 = 200. Because we recover General 
Relativity on large scales (e — > for k — > 0) all curves 
converge to the ACDM threshold at large mass and show 
increasingly large deviations from GR at smaller mass. 
The asymptotic value is 5 C ~ 1.59 rather than 1.67 as we 



define 5 C as T~ 1 (20Q) instead of J r g " 1 (oo), that is, by a 
nonlinear density contrast of 200 rather than by the full 
collapse to the center, as in [65]. 

Similar trends were obtained in [2l| . using a simpli- 
fied dynamics described by an effective Newton constant 
that depends on the "environment" density, which al- 
lowed them to include screening effects. Thus, because 
the latter are more important for large mass they ob- 
tained a mass-dependent threshold 5 C that decreases at 
small mass and converges to the GR value at large mass. 
We can see in Fig. [5] that even without such screening 
effects, a dependence on mass is already present because 
of the dependence on wavenumber of e(/c,a). Since both 
effects show similar trends, including them both would 
give a steeper dependence on mass than in Fig. [5J Nev- 
ertheless, it is interesting to also investigate both mecha- 
nisms separately, as their relative amplitude depends on 
the details of the modified-gravity model. 

We can see that the approximation (|102j) somewhat un- 
derestimates the departure from the GR result. This can 
be understood from the fact that the dynamics steepens 
the density profile, which amplifies the right hand side in 
Eg . (l9"6")l . Nevertheless, the approximation (|102l) . which 
is much easier to compute, gives a reasonable estimate 
of the modified-gravity effect. Because inner shells have 
already collapsed when the shell at mass M reaches the 
nonlinear threshold S x = 200, we should include virial- 
ization effects which smooth out the inner density profile. 
Therefore, the difference seen in Fig.[5]should actually be 
somewhat overestimated. Moreover, for smaller nonlin- 
ear density contrast S x the relative deviation decreases, 
because the ansatz (llOOp is exact at linear order (for our 
initial conditions). Thus, for practical estimates the ap- 
proximation (|102p should be sufficient, at least in a first 
step. 

VI. DENSITY CONTRAST PROBABILITY IN 
THE QUASILINEAR REGIME 

Following (6f| l66j . we can use the spherical collapse 
dynamics described in Sect. [V]to derive the probability 
distribution of the matter density contrast in the quasi- 
linear regime. 

To compute the probability distribution, V(S X ), of the 
nonlinear density contrast within a sphere of comoving 
radius x, it is convenient to introduce the cumulant gen- 
erating function 

e -v(y)/^ (A) = / e -v*./°2 C A)\ (106) 

dS x e -^/<A) V(S X ). (107) 

This determines the distribution V{5 X ) through the in- 
verse Laplace transform 

/+ioo 1 
V — e bfe-v(y)]/^ (A ). (108) 
-ice 2 ™ X (A) 
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In Egs. (jl07|) - (|108|) we rescaled the cumulant generating 
function by a factor & x t A \ so that it has a finite limit in 
the quasilinear regime, QWa) — > 0, for the case of Gaus- 
sian initial fluctuations [381 ] . In particular, its expansion 
at y = reads 



v{y) 



E 



(-y) n (^)c 



2(n-l) 



(109) 



'x(A) 



The average (I106P can be written as the path-integral 

e -v<v)/ol w = (detCr 1 . (A) ) 1/2 hh(A) e~ s ^^ 



where ^Sl5l(A) 1S ^ e mverse matrix of the two-point 
correlation of the reference linear density field and the 
action S reads as 



S[S L (A)} = V$x[$L(A)] 



'x(A) 



?L(A) 



' C S L S L {A) ■ 5 L(A) 

(HI) 

Here ^[^(A)] is the nonlinear functional which assigns 
to the initial condition, defined by the reference linear 
density field <5l(a)(x'), the nonlinear density contrast 5 X 
within the sphere of radius x. 

As in Sect. IV C| we choose to define the initial con- 
ditions through the reference ACDM linear field $l(a)- 
We could also write all expressions above in terms of the 
actual linear field 6l, its correlation Cs L $ L , and the vari- 
ance u x . Here we prefer the formulation (jllOp because it 
clearly separates the initial conditions from the modificd- 
gravity effects. Thus, in the action (jlll[) all modificd- 
gravity effects are enclosed in the functional <$ x [£.l(a)]j 
which describes the gravitational dynamics, whereas if we 
express the initial conditions in terms of the e-dependent 
linear field 5l these modified gravity effects would appear 
in all terms of the action. Of course, we adopt this for- 
mulation because we wish to compare with this ACDM 
reference several models that only show small deviations. 

The action S does not depend on the normalization of 



the linear power spectrum since both trLj 



and C, 



SlSl(A) 



are proportional to Pl(A)- Then, in the quasilinear limit, 
fx(A) "~ * 0, the path integral (|110l) is dominated by the 
minimum of the action 66], 



r x (A) : ip{y) -> 



)(x') 



S[8 L {A)\- 



(112) 



Using the spherical symmetry of the top-hat window W 
that defines the spherical average S x , one obtains a spher- 
ical saddle-point [66J. In General Relativity its linear 
radial profile is given by Eq. (|104l) . where q is the La- 
grangian radius that corresponds to the Eulerian radius 



x, 



q 3 = (l + S x )x 3 



(113) 



Then, the amplitude £>l<?(A) of the saddle-point (|104l) . 
which also sets the scale q through Eq. (|113p , is given 
by the spherical-collapse mapping, 



This derivation agrees with the results that can be ob- 
tained from a perturbative computation of the cumulants 
(S x ) c at leading order and a resummation of the series 
(|109[) (gH . It also extends these results to the case where 
the series (|109l) has a zero radius of convergence, which 
occurs when V(5 X ) decreases more slowly than a simple 
exponential at large densities (6(| 

A nice feature of this derivation is that it bypasses 
the computation of the cumulants {5 X ) C through the ker- 
nels of Eq. (f52|) , as all spherically-averaged quantities 
are given by the spherical-dynamics mapping .F(<5l(a)) 
(which includes terms at all orders by expanding over 
&l(A))- However, the problem is more complex in our 
case because of the dependence of e(k, a) on wavenumber. 
Indeed, this means that the nonlinear density contrast 6 X 
at radius x does not depend on the linear density con- 
trast <5l 9 (a) at the Lagrangian radius q, associated with 
the same mass M only. Indeed, as discussed in Sect. |V] 
the spherical dynamics (|9"6"| depends on the full shape of 
the initial perturbation. Taking into account this modi- 
fication changes the profile <$l(a)( x ') of the minimum of 
the action S[Si,r A -\] m Eq. (|112p , because the functional 
<y#i(A)(x')] is no longer of the form 5 X = J"(<Wa))- 

To simplify the analysis we neglect this change of the 
profile of the saddle-point. This is actually valid to 
first order over e. Indeed, let us write the action S 
as S = So + eSi, where So is the usual ACDM ac- 
tion (where e = 0), and S± is the modification due to 
a nonzero e(k, a) kernel, where we factored out a nor- 
malization parameter e that scales as e. Because of 
this new term eSi, the saddle-point £>l(A) is changed to 
<Sl(A) = ^LQ (A) + e^Li(A), where 5 L0 (A) is the GR saddle- 
point (|104p . Then, the generating function is changed to 
<P(y) -> S [S L0 {A) + e^Li(A)] + e<Si[£io(A) + eS L1(A) }. Be- 
cause 5lo(a) is a saddle-point of the action So, we have 
So[ho(A)+£dLi(A)} = S , o[<5lo(A)]+C(£ 2 ), that is, S [S LW ] 
is only modified by terms of order e 2 . Because of the pref- 
actor e we also have eS , i[(5 L0 (A) +^li(A)] = eS r i[5 i0 (A)] + 
0(e 2 ). Therefore, S[S L(A) ] = S[S L0(A) } + 0{e 2 ) and we 
can neglect the change of the saddle-point up to first or- 
der over e. In fact, we do better than this because we 
only neglect the change of the radial profile but we keep 
track of the dependence on e of the amplitude &Lq(A) °f 
the saddle-point. 



On the other hand, if we use the approximation (|102p 
instead of Eq. ([M|) . the functional <y<5L(A)( x ')] is again 
of the form S x = J~ q {8L q (A)) an d the saddle-point profile 
(|104p becomes exact within this approximation. 



5 X = J 7 (5 i9 (A))- 



(114) 



In both cases, whether we use the approximation (|102p 
or the exact equation ([TO]) , the function T q now also de- 
pends on the scale q, in contrast to the usual Newtonian 
case. 

Then, from this spherical-collapse mapping ^(^l^a))) 
described in Sect. V C[ we obtain the generating function 
(p(y) as follows [65, 6f|. Substituting the profile (|104p into 
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the minimum (|115l) also writes as 



FIG. 9: Probability distribution of the matter density con- 
trast within spherical cells of radius 5h~ Mpc at z — (all 
curves almost fall on each other). 
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FIG. 10: Relative deviation from General Relativity of the 
probability distribution V(5 X ), at redshift z = for a radius 
x = 5/i _1 Mpc. For each (n, mo) model the deviation from 
GR is positive at low and high densities and negative around 
S ~ 0. The solid and dotted lines are the exact results from 
Eq.® for (n,m ) = (1,0.1) and (0,0.1). The closest dashed 
line of the same color is the result from the approximation 
(11021) . for the same value of (n, mo). 



Eq. (jllip and using Eq. (j!14j) the minimum (|112l) reads as 

(115) 



(f(y) = min 



9(A) 



2 a z ~ LqW 



Defining the function t(£) through the parametric system 

MM, 



C = 5 X = J"g(^i,g(A)) and t = S Lq{A y 



C (A) 



Tq(A) 



(116) 



f(y) — min 



r(C) 2 



(117) 



This corresponds to the implicit equations (Legendre 
transform) 

dr r 2 
V = -r-^ and ip = y( + —. (118) 

Finally, this gives the probability distribution V(S X ) 
through Eq. (fT08l) . The probability distribution V(S X ) 
depends on the spherical-collapse dynamics and on the 
shape of the initial power spectrum P L i A){k) 1 through 
the ratio cr a; (A)/ a 'q(A) in the second Eq. (|116j) . This sec- 
ond effect, sometimes called a "smoothing effect" jsij], is 
due to the collapse (or expansion) of the mass shell M 
from the Lagrangian scale q to the Eulerian scale x. This 
mixes scales and implies that the distribution V(S X ) at 
scale x is sensitive to the initial power over all scales. In 
our modified-gravity second dependence on the 

shape of the linear power spectrum appears through the 
mapping T q itself, because of the e-dependent terms in 
Eqs.dM]) and (fTU2|) . 

We show in Fig.|9]the probability distribution ViS^) at 
redshift z — and radius x = 5/i _1 Mpc. Here we use the 
exact dynamics (|96[) but using the approximation (|102[) 
gives very close results that would not be distinguished in 
this figure. We recover the usual asymmetric shape due 
to nonlinear gravitational clustering, which builds an ex- 
tended high-density tail and shifts the peak of the distri- 
bution towards low densities before a sharp low density 
cutoff at 8 X — > — 1 + (on small scales, most of the mat- 
ter lies in overdensities but most of the volume lies in 
underdense regions). 

Since it is difficult to distinguish different curves on this 
figure we plot the relative deviation from GR in Fig. [TOJ 
for the two models where it is the largest. (The two other 
cases would fall below the range plotted in the figure for 
the most part.) We plot our results using either the ex- 
act equation or the approximation (|102[) . We can 
see that both curves are very close. Indeed, as explained 
in Sect. IV CI for smaller density fluctuations the ansatz 
(jlOOP becomes more accurate as it is exact to linear order 
and the profile has not had time to be strongly modified 
by the dynamics (moreover, the collapse is not very sen- 
sitive to the exact shape of the profile). 

As we consider models with a positive value of e, which 
leads to an effective amplification of gravity, it is easier to 
build large nonlinear density fluctuations. This was also 
apparent in Fig. [S]for the specific case of S x — 200. For 
Gaussian initial conditions the tails of the probability dis- 
tribution ViSx) are of the form V^x) ~ e~ 



'Lq(A) 



/(2<A)) 



where <^l 9 (a) = J~ q ~ 1 (S x ), and the lower value of |<5l 9 (a)| 
that is needed to reach a given \S X \ yields a slower decay 
of the rare-event tails. This is why we recover a positive 
deviation from GR (i.e., a higher probability P) at both 
very low and very high densities in Fig. 1101 Of course, 
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since probability distributions are always normalized to 
unity this implies that the relative deviation shows a 
change of sign and that the probability distribution V(5 X ) 
obtained in these models is smaller than the ACDM one 
for moderate densities. This explains the behaviours seen 
in Fig. EO] 

These features are in qualitative agreement with the 
results obtained in numerical simulations of various mod- 
ified gravity models [H, ISSl ; which also find that an effec- 
tive amplification of gravity generically leads to more nu- 
merous very low density and high density regions, while 
shifting the peak of the probability distribution towards 
lower densities. 

The relative deviation from GR does not necessarily 
grow to unity at high densities (and may even decline). 
This is due to the fact that high densities at a given Eu- 
lerian radius x correspond to large masses, hence to large 
Lagrangian (i.e. initial) radius q. Then, because we re- 
cover General Relativity on large scales the linear thresh- 
old <5l 9 (A) = J'q 1 {Sx) converges to the one obtained in 
the ACDM cosmology, as in Fig. [SJ Therefore, depend- 
ing on the rate of convergence towards General Relativity 
on large scales (as compared with the increasingly high 
sensitivity of the rare tail) the large-density tail may or 
may not converge back to the GR prediction. In mod- 
ified gravity scenarios with a screening mechanism that 
implies convergence to GR in high-density environments, 
such as the chameleon mechanism, the high-density tail 
is expected to show a faster convergence back to the GR 
prediction. 

These effects do not appear at very low densities, which 
correspond to increasingly small mass M and Lagrangian 
radius q, where the modifications from General Relativ- 
ity do not vanish within our framework. In this limit, 
the relative deviation of V(S X ) from the ACDM reference 
can grow up to unity. However, this appears far in the 
low-density tail, which is characterized by a very sharp 
cutoff, and this may not be a very efficient tool to probe 
modified-gravity effects. 



VII. HALO MASS FUNCTION 

The computation of the probability distribution V(S X ) 
was described in the previous section for the quasilinear 
regime, ct^a) — > 0. However, this result is more general 
and actually applies to rare events, where the path inte- 
gral (|110[) is peaked around the minimum of the action 
S. In the quasilinear limit any finite nonzero density con- 
trast 5 X becomes a rare event, which is why Eq. (|117[ ) de- 
termines the full probability distribution in this regime. 
For arbitrary values of a x , Eq. (|117l) applies to rare events, 
that is, to the tails of the probability distribution V{8 X ) 
[sU (this again allows one to recover the results obtained 
from a perturbative analysis (63[). However, for large 
overdensities shell crossing appears at some stage (typi- 
cally for S x > 200), after which Eq. (|117[) no longer holds 
[64ll65l|. Nevertheless, for lower densities one obtains the 
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FIG. 11: Halo mass function at redshift z — 0. 
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FIG. 12: Relative deviation from ACDM of the halo mass 
function at redshift z — 0. 



asymptotic behaviour V{8 X ) ~ e _ **«C A >/( 2o «CA)). This also 
determines the large-mass tail of the halo mass function 
n(M)dM/M, where we define halos as spherical objects 
with a fixed density contrast threshold 6 — 200, 



M 



ln[n(M)] - - 



2a (A) (M)2 



with 



6 L(A) (M)=F- 1 (d), 



(119) 



(120) 



where at^(M) = <r g (A) with M = p m 4irq 3 /3. 



As in |6lL l65| . a simple approximation for the mass 
function that satisfies the large-mass asymptote (IL L9[) can 
be obtained using the Press & Schechter scaling variable 

v [z3, 



v M M v 



(121) 
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with 



V(200) 
'(A)(M) ' 



(122) 



where we choose to define halos by the nonlinear density 
threshold 8 — 200. The scaling lunction f(v) is obtained 
from a fit to ACDM numerical simulations that satisfies 
the exponential tail f(v) ~ e~ v ' 2 [65| 



f{v) = 0.502 [(0.6z/ 



2.5 



(0.62i/) 



0.51 



- V 2 /2 



(123) 



This ensures that the halo mass function is always nor- 
malized to unity and obeys the large-mass tail (I119[) , for 
any spherical-collapse mapping T q . The only change 
from the ACDM cosmology is that the linear thresh- 
old J" (J _1 (200) in Eq. (|122l) now depends on the mass M 
through the scale q(M). The approximation (|123p only 
ensures that the large-mass tail is correct, but it may 
happen that the low-mass power-law tail should depend 
on e. An analysis of such effects would require numer- 
ical simulations because analytical methods cannot pre- 
dict the low-mass tail of the halo mass function (which is 
sensitive to mergers and non-local effects). Nevertheless, 
we can expect modifications for moderate masses to be 
less important and partly taken into account through the 
normalization constraint of the mass function. 

As compared with the excursion set approach pre- 
sented in [2ll . I711 1721 ] , we do not include screening effects 
but we take into account the dependence on wavenum- 
ber of the modified-gravity kernel e(k, a). As explained 
in Sect. El this leads to a mass-dependent linear thresh- 
old Sl(M) whence to deviations from the ACDM mass 
function that will depend on mass. 

We show the halo mass function in Fig. [TTJ and its rel- 
ative deviation from the ACDM mass function in Figfl^l 
Here we use the approximation (|102l) for the mapping 
3~q{&Lq(h)) but we checked that using Eq. (|96)l yields close 
results. For the models that we consider here the mass 
functions are very close to each other and relative devia- 
tions are on the order of 10% or less. In agreement with 
the behaviour of the probability distribution V(S X ) dis- 
cussed in the previous section, a positive e(fc, a) leads to 
more numerous high density fluctuations and to a larger 
number of massive collapsed halos. This explains why the 
ratio to the ACDM mass function is greater than unity 
for v > 1, which corresponds to rare halos. Again, this 
relative deviation grows for lower n and smaller mo- 

The same trends appear in numerical simulations of 
similar modified gravity scenarios [H, [68|, [73|, , with an 
increase of the large-mass tail for models with an effective 
amplification of gravity. We show our results for f(R) 



models with \f RQ \ = 10~ 4 , 10~ 5 , lO" 6 , as in [2J, H UA 
in Appendix [A"l 

On the mass scales shown in FigfT2l the ratio keeps 
growing at high masses for mo = 0.1 while it decreases 
for mo = 1- As in the high-density tail shown in Fig. [TOl 
this is due to two competiting effects: i) the exponential 
tail (|119l) of the halo mass function amplifies the sen- 
sitivity to modified-gravity effects at large masses, but 



ii) these deviations from General Relativity decrease at 
large scale whence at large mass (e(k, a) — > for k — > 0), 
as seen in Fig. [3] Then, depending on the relative im- 
portance of both effects, the ratio of the mass function 
to its ACDM reference may or may not grow with mass 
on the scales that are considered. As expected, a lower 
parameter too (which implies a modification of gravity 
up to larger scales, k ~ mo and q ~ 1/too, see Eq.fJTT])) 
yields a slower convergence to General Relativity at high 
mass, whence a larger weight to the first effect i) above. 
This explains why on the mass scales shown in FigfT2lthc 
ratio keeps growing at high masses for too = 0.1 while it 
decreases for too = 1. 



VIII. FROM LINEAR TO HIGHLY 
NONLINEAR SCALES 

Following [33, IH| > we can combine the perturbative re- 
sults of Sect. II VI with the halo mass function of Sect. IVIII 
to obtain the matter density power spectrum and bis- 
pectrum from linear to highly nonlinear scales. As in 
the usual halo model [ill, we write the nonlinear power 
spectrum as the sum of "two-halo" and "one-halo" terms, 



P(k) = P 2H (k)+P m (k), 



(124) 



where P2H is the contribution from pairs of particles that 
are located in two different halos and Pm is the contri- 
bution from pairs located in the same halo. As explained 
in |37| , P2H contains the perturbative contribution to the 
power spectrum and we write 



P2n(k) =P 2H (27r/fc)P pert (fc), 



(125) 



where P2h(<z) is the fraction of pairs, with initial (i.e. La- 
grangian) separation q, that belong to two distinct halos, 
and Ppert(fc) is the power spectrum obtained by pertur- 
bation theory. It is not possible to use the standard one- 
loop prediction, unless one adds a high-fc cutoff, because 
it grows too fast at high k and leads to unphysical re- 
sults at high k for the sum (|124[) . Here we consider the 
one-loop prediction P pC rt(fc) given by the resummation 
with Eqs. (l6"6"|) - (f6"7} . Indeed, at this order it yields 



Ppert(fc) ~ Ph{k) at high k [33j, so that the two-halo 
term is subdominant with respect to the one-halo term 
and one obtains a good match to numerical simulations 
[37l |45|. Next, the one- halo contribution, which is fully 
nonperturbative, reads [13] 

PlH(fc) = I" — /M^"T3 (uM{kf - W(kqf) , 
Jo v P m (27r) J V / 

(126) 

where W is the Fourier transform of the 3D top-hat, de- 
fined in Eq. (1901) , and % is the normalized Fourier trans- 
form of the density profile pm(x) of halos of mass M, 

u M (k) = JL y dx e- ik ' x p M {x). (127) 
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We use the usual "NFW" halo profile (76| , with the mass- 
concentration relation from [37j. Therefore, we do not 
take into account the effects of the modified gravity on 
the shape of the profiles of the dark matter halos. Our 
one-halo term Pm only depends on e(fc, a) through the 
change of the halo mass function described in Sect. IVIII 
The counterterm W 2 in Eq. (|126p ensures that the one- 
halo contribution decays as Pm(k) oc k 2 at low k, so 
that the total power (I124[) converges to the linear power 
on large scales. This follows from the conservation of 
matter and the fact that halo formation corresponds to 
a small-scale redistribution of matter [37], [z3 [sH ■ 

In a similar fashion, the matter density bispectrum can 
be written as the sum of three-halo, two-halo, and one- 
halo terms, 

B = B m + B 2li + B m: (128) 

with m, 

B 3 n(h, k 2 , k 3 ) = B peit (ki, fc 2 , k 3 ), (129) 
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FIG. 13: Logarithmic power, A 2 (k) = 47rfc 3 P(fc), at z = for 
four (n, mo) models. In each case we plot the linear power 
(dashed line) and the nonlinear power (solid line). 



B 2 n(ki,k 2 ,k 3 ) = P L (h) J 



Av M 



fiy) 



x Y[ (uuikj) - Wfaq)) + 2 eye, (130) 

3=2 ' 



Bm(ki,k 2 ,k 3 ) 

3 



dv t, \ ( M 



II - W{k 3 q)^ + 2 eye, (131) 



3=1 



Again, the counterterms W in Eqs. (|130p and (1131|) en- 
sure that the two-halo and one-halo contributions decay 
on large scales so that the bispectrum converges to the 
perturbative prediction -B pcr t- As found in [45| and con- 
trary to the situation encountered for the power spec- 
trum, the standard one-loop perturbation theory predic- 
tion for -Bp Cr t is well-behaved at high k (i.e., it is signif- 
icantly smaller than the one-halo contribution) and it is 
more accurate than the resummation schemes that have 
already been studied. Therefore, we only consider the 
standard perturbative approach for the three-halo contri- 
bution (|129p . More precisely, we use the exact tree- level 
result (|T9")) and the approximate one-loop correction ([50)1 
by setting B peTt = B tICC + B lloop . 

While Eq. (I128[) yields a reasonably good match to nu- 
merical simulations (~ 10%) over all scales for the bis- 
pectrum [45j], Eq. (|124p significantly underestimates the 
power spectrum on the transition scales (by ~ 20 — 30%), 
even though it gives a good accuracy on larger scales 
(~ 1% below k ~ 0.3/i Mpc -1 at z = 1) and smaller 
scales (~ 10% above k ~ 5h Mpc" 1 at z = 1). Fol- 
lowing I l , we consider a simple power- law interpolation 
^tangbetween large and small scales, 

J Ptan g (fc) = P 2 n+m(k) for k < fc_ and k > k' + (132) 
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FIG. 14: Equilateral bispectrum B Bq (k) = B(k, k, k), at 2 = 
for four (n, mo) models. Fhe bispectrum is multiplied by a 
factor k 3 in this plot to decrease the range spanned by the 
vertical axis and to make the figure easier to read. In each 
case we plot the tree-level bispectrum (dashed line) and the 
full nonlinear bispectrum (solid line). 



and 



^tang(fc) is a power law within fc_ < k < k' + . (133) 

The transition range [fc_, k' + ] is automatically determined 
from the shape of P2H+in(k) and B(k,k,k) and it de- 
pends on the shape of the linear power spectrum and 
on redshift. This improves the agreement with numeri- 
cal simulations in the ACDM cosmology [45| while keep- 
ing the perturbative and 1-halo behaviours on large and 
small scales. 

We show in Figs. 1131 and H~4l the matter density power 
spectrum and bispectrum that we obtain at redshift 
z = 0, from linear to highly nonlinear scales. The var- 
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FIG. 15: Relative deviation from ACDM of the power spec- 
trum obtained in four models at redshift z = 0. In each case, 
we plot both the relative deviation of the linear power (dashed 
line) and of the nonlinear power (solid line). From left to right 
we consider the models (n, mo) = (0,0.1), (1,0.1), (0, 1), and 
(1,1). 
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FIG. 16: Relative deviation from ACDM of the bispectrum 
obtained in four models at redshift z — 0. In each case, 
we plot both the relative deviation of the tree-level bispec- 
trum (dashed line) and of the nonlinear bispectrum (solid 
line). From left to right we consider the models (n, mo) = 
(0,0.1), (1,0.1), (0,1), and (1,1). 



converge to the GR behaviour. Then, we expect that 
our modclization provides a similar accuracy to the one 
found in ACDM cosmology by comparison with numeri- 
cal simulations [37l. l45j. As in Sect. IIVD1 we clearly see 
that nonlinear gravitational clustering amplifies both the 
power spectrum and bispectrum at hig h k but damps the 
baryon acoustic oscillations. As in |37l |45| . our approach 
allows us to describe the power spectrum and bispectrum 
from large linear scales down to small highly nonlinear 
scales. 

We show in Figs. [15] and [16] the relative deviations 
from the ACDM reference of the power spectrum and 
of the equilateral bispectrum. In the weakly nonlinear 
regime the relative deviations grow with k, following the 
behaviour of e{k,a). In agreement with the discussions 
above, they reach a maximum on transition scales, start- 
ing to deviate from the ACDM growth for k ~ too, and 
then slowly declining on highly nonlinear scales. On these 
nonlinear scales, the relative deviations at the level of the 
linear or tree-order contributions are no longer a good es- 
timate of the actual signal and greatly overestimate the 
effects of modified gravity. Since the theoretical accuracy 
is greater on weakly nonlinear scales (which can be an- 
alyzed by systematic perturbative approaches) than on 
highly nonlinear scales (which require phenomenological 
ingredients such as halo profiles), these behaviours sug- 
gest that it is more efficient to focus on weakly nonlinear 
scales to probe such modifications of gravity. 

It is is also worth emphasizing that the deviations 
from ACDM which we have calculated with the steep- 
est descent resummation method together with the halo 
model show the same trends as the N-body results 
25] obtained for models with n — 1 and |/ro| = 
10~ 4 , 10~ 5 , 10~ 6 . Indeed, numerical results show that 
the deviation from A-CDM reaches a peak at weakly non- 
linear scales before decreasing on highly non-linear scales. 
Simple fitting procedures designed for ACDM cosmology 
[78j have been shown not to provide good results and to 
miss this high-fc behavior [25|]. This shows the advan- 
tage of approaches like ours that are closer to physical 
modeling. Even though they may be less accurate than a 
specific fitting formula , their behaviour as cosmological 
parameters and scenarios are modified is more reliable. 



IX. CONCLUSION 



ious curves are very close and we can see that at high 
k the deviations are actually damped by nonlinear ef- 
fects. Within our framework, this is because we ne- 
glected any impact of modified gravity on the halo profile 
(|127[) and the only influence of modified gravity appears 
through the halo mass function n(M). This may not 
be such a bad approximation because in more realistic 
models modifications to gravity vanish on small scales 
(e.g., through chameleon or Vainshtein mechanisms) so 
that the density profiles of small halos are expected to 



We have considered the dynamics of structure forma- 
tion in modified gravity models analytically. To do so, 
we have used a steepest descent technique for the gen- 
erating functional of density and velocity perturbations 
as well as the spherical collapse dynamics. The models 
we have considered correspond to screened modifications 
of gravity due to a scalar field. In numerical examples 
we have focused on models defined by a power law mass 
function and a constant coupling to matter, which coin- 
cide with f(R) models in the large curvature limit and in 
the matter era, although the techniques developed here 
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are general. The results we have presented comprise the 
power spectrum, the bispectrum, the probability distri- 
bution of the density contrast, and the large-mass tail of 
the halo mass function. Modified gravity has interest- 
ing features astrophysically when the ratio of the mass 
of the scalar field over the Hubble rate now too / Hq is of 
order 10 3 . In this case, deviations can be substantial and 
larger than a few percent. In this paper, we do not at- 
tempt to give precise predictions, we are more interested 
in indications that can be obtained relatively fast using 
our analytical tools without the need for large N-body 
simulations. 

After a description of the linear growing and decay- 
ing modes, which become fc-dependent in these modificd- 
gravity scenarios, we have obtained the associated linear 
growth rate f(k,z). For the realistic parameters (n, mo) 
studied here measuring its deviation from the General 
Relativity prediction remains challenging, but future sur- 
veys such as Euclid should give a clear signal for the most 
favorable cases (e.g., (n, mo) = (0,0.1)). 

Next, we have described how higher-order perturba- 
tive contributions can be computed in the weakly non- 
linear regime. The dependence on wavenumber of the 
linear modes makes numerical implementations of these 
perturbative schemes significantly more complex than in 
the usual General Relativity case, because time- and 
scale-dependences no longer factor out. We have pre- 
sented the generalization of the "standard" perturbative 
approach as well as a "steepest descent" approach that 
performs partial resummations of higher-order diagrams. 
The path-integral formalism that underlies this second 
method also provides an efficient route to recover the 
standard perturbative approach and avoids the need to 
compute the n— point kernels F£. We find that for re- 
alistic modificd-gravity scenarios, such as the ones in- 
vestigated here, the deviations of the power spectrum 
from General Relativity on BAO scales are quite modest 
(typically less than 6%) and below the accuracy of the 
standard perturbative approach at one-loop order. This 
means that one must use more accurate schemes, such as 
the one-loop steepest-descent approach presented here, or 
possibly include higher-order terms within the standard 
approach (but its convergence is not very well behaved). 

For the bispectrum we find that nonperturbative con- 
tributions (associated with one-halo and two-halo terms) 
cannot be neglected on the weakly nonlinear scales where 
the deviations from General Relativity can be detected. 
This suggests that for practical purposes the power spec- 
trum is a more reliable probe of such modified-gravity 
effects, because its deviations from the GR predictions 
are larger than for the bispectrum in the perturbative 
regime, where rigorous and systematic approaches can 
be developed. 

To go beyond these low-order perturbative approaches, 
we have described the dynamics of spherical density fluc- 
tuations, which can be exactly solved before shell cross- 
ing. Again, modifications to gravity make the analysis 
significantly more complex, because the motions of differ- 



ent shells no longer decouple, even before any shell cross- 
ing. This means that one must solve the evolution with 
time of the full density profile. Nevertheless, we have in- 
troduced a simple approximation for typical profiles that 
allows to decouple the motion of the mass shell of inter- 
est. We find this provides a reasonable approximation 
to the exact dynamics (but slightly underestimates the 
effects of modified gravity). This analysis provides the 
characteristic dependence on mass of the critical linear 
density threshold 5 C (M) associated with a given nonlin- 
ear threshold (such as S — 200). In the cases studies here, 
where the function e(fc, a) is positive and corresponds to a 
time- and scale-dependent effective amplification of grav- 
ity, this threshold S C (M) decreases at low mass (because 
this amplification is larger on smaller scales) and con- 
verges to the constant GR prediction at large mass (be- 
cause we recover General Relativity on large scales). 

In contrast to some previous works, this dependence on 
mass does not arise from screening effects (that depend 
on mass through the depth of the gravitational potential, 
which triggers the screening mechanism) but from the In- 
dependence of the modified-gravity kernel e(k, a). 

This also allows us to obtain the probability distri- 
bution, V(S X ), of the nonlinear density contrast within 
spherical cells, in the weakly nonlinear regime. Because 
of this effective amplification of gravity, the tails of V(5 X ) 
grow with respect to the General Relativity prediction 
(and by conservation of the probability normalization to 
unity V{8 X ) decreases for moderate density fluctuations). 
This growth is smaller and the relative ratio to GR does 
not necessarily goes to infinity in the large-density tail, 
as opposed to the low-density tail, because on large scale 
the dynamics converges to General Relativity. 

The same effect amplifies the large-mass tail of the halo 
mass function. Again, the ratio to the GR prediction may 
increase or decrease with mass in the range of interest 
depending on how fast modifications to gravity vanish 
on large scales. 

Finally, combining perturbative approaches with halo 
models, we have computed a simple estimate of the power 
spectrum and bispectrum from linear to highly nonlinear 
scales. Within this modclisation, we find that the rela- 
tive deviation from General Relativity is the largest on 
the transition scales between the linear and the highly 
nonlinear regimes, for both the power spectrum and bis- 
pectrum. Since nonlinear scales are difficult to predict 
with a high accuracy (because of the complex nonper- 
turbative dynamics associated with shell crossings and 
because one should include baryon and galaxy formation 
effects), this suggests that weakly nonlinear scales, in par- 
ticular in the perturbative regime, are the best probes of 
these modified-gravity models. 

As a summary, our new results can be listed as follows: 

- a comparison of the accuracy of one- loop perturbative 
expansions (by using two such schemes and by estimat- 
ing non-perturbative one-halo contributions) with realis- 
tic deviations from GR, for the matter power spectrum 
and the bispectrum. 
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- an analysis of the spherical collapse that includes 
shell-coupling and the scale-dependence of the modificd- 
gravity kernel e(fc, a). 

- the dependence on mass, due to the scale-dependence 
of e (and not to screening effects), of the deviation from 
GR on the halo mass function. 

- an analytical model for the probability distribution 
V(S X ) in the rare-event regime. 

- a combination of one-loop perturbative expansions 
with halo models for the matter power spectrum and the 
bispectrum up to highly non-linear scales. 

Our methods call for improvements to reach the needs 
of precision cosmology. Indeed we have neglected, for 
ease of treatment and as a first step, two major effects. 
The first one consists in including non-linearities in the 
scalar field sector of the models. Here the scalar field dy- 
namics are only linear and non-linear effects in both the 
potential and the coupling to matter ought to be consid- 
ered. Technically, this can be done at the one loop level 
by self-consistently modifying the Euler equation with 
non-linear terms coming from the scalar field interaction 
with matter particles. A second ingredient we have not 
considered so far is the screening of the scalar force in 
dense environments. This will modify the spherical col- 
lapse of over densities and therefore the halo statistics. 
Eventually this will have an impact on the growth of 
non-linear structures. As a result, the effects described 
in this paper can only be taken as indications on quasi- 
linear scales. Work on all these aspects is in progress. We 
also intend to carry out a comparison of our analytical 
results with the N-body simulations which use the same 
mass and coupling parameterisation of modified gravity. 
Doing so, and for a greater variety of models including 
dilatons and symmetrons, we hope to validate our an- 
alytical approach which could then be used for models 
that will appear in the future and be analysed without 
the need for large N-body simulations. 



Appendix A: The case of f(R) models 

We consider in this appendix the case of f(R) models 
which have also been studied through numerical simula- 
tions, with a power-law form as in Eq. (II6p . The mass of 
the scalar field evolves with time as [23 



m(a) = mo 



Ono(l + z) 3 + 4«ao 



4«A0 



(n+2)/2 



(Al) 



where mo is given by Eq. (|17[) . This gives the approximate 
scaling (p~5|) at high redshift but for accurate computa- 
tions it is necessary to use the more precise expression 

CEO). 

To compare with the numerical results of [U, [25|, 
Izl III Izll we adopt the same WMAP3 ACDM 
reference model |80j, with cosmological parameters 
(Q m ,n h ,h,a 8 ,n B ) ' = (0.24,0.04181,0.73,0.76,0.958). 
We focus on the case n = 1, with the amplitudes 
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FIG. 17: Relative deviation from ACDM of the halo mass 
function at redshift 2 = 0, for n = 1 and |/fl | = 10 -4 , 10 -5 , 
and 10 -6 , from top to bottom. 
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FIG. 18: Relative deviation from ACDM of the power spec- 
trum at redshift 2 = 0, for n = 1 and |/a | = 10~ 4 ,10~ 5 , 
and 10 -6 . In each case, we plot both the relative deviation 
of the linear power (dashed line) and of the nonlinear power 
(solid line). The points are the results of the "no-chameleon 
simulations" from 12511. 



|/ flo | = 10" 4 , 10~ 5 , and 10~ 6 . We show in Figs. [13 M 
and [T9l the relative deviations from the ACDM refer- 
ence of the halo mass function, the matter power spec- 
trum, and the bispectrum. Our results are similar to the 
ones obtained in Figs. HS1 and [TBI in the main text, 
for our power-law models parameterized by (n,mo). We 
can check that our results also show a reasonable agree- 
ment with the "no-chameleon" numerical simulations of 
[H, [zl, [ll] for the halo mass function and the power 
spectrum, although we may overestimate the large-mass 
tail for M > lO 15 /i~ 1 Af . The almost straight lines on 
transition scales in Fig. [15] correspond to the interpo- 
lation (|133|) and should not be considered as an accu- 
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FIG. 19: Relative deviation from ACDM of the bispectrum 
at redshift z = 0, for n — 1 and |/i? | = 10 _4 ,10 -5 , and 
10 -6 . In each case, we plot both the relative deviation of 
the tree-level bispectrum (dashed line) and of the nonlinear 
bispectrum (solid line). 

rate prediction. However, they correctly reproduce the 
saturation of the relative deviation and the transition 
towards the highly nonlinear regime (dominated by the 



one-halo contribution) where the relative deviation de- 
clines (within our framework, where we neglect modifi- 
cations of halo profiles) . The same behaviour is found in 
numerical simulations, as can be seen in figure 18 where 
we compare our results to the no-chameleon simulation 
of 25], with a reasonably good quantitative match. It is 
interesting to note that using simple fitting procedures 
designed for ACDM cosmology, such as the halo-fit from 
78], has been shown not to provide good results and 
to miss this high-fc behavior [25|]. This is not fully sur- 
prising, since such fitting formulae were not designed for 
these scenarios. This shows the advantage of using ap- 
proaches such as the one presented in this paper that are 
closer to physical modeling (using both systematic per- 
turbative expansions and phenomenological halo mod- 
els) . Even though they may be less accurate than a spe- 
cific fitting formula for the class of models the latter was 
built from, their behaviour as cosmological parameters 
and scenarios are modified is more reliable. 
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